March 12, 2015

**[This is a guest post by Eric-Jan Wagenmakers and Quentin Gronau introducing the RGraphCompendium. Click here to see the full compendium!]**

Every data analyst knows that a good graph is worth a thousand words, and perhaps a hundred tables. But how should one create a good, clean graph? In R, this task is anything but easy. Many users find it almost impossible to resist the siren song of adding grid lines, including grey backgrounds, using elaborate color schemes, and applying default font sizes that makes the text much too small in relation to the graphical elements. As a result, many R graphs are an aesthetic disaster; they are difficult to parse and unfit for publication.

In constrast, a good graph obeys the golden rule: “create graphs unto others as you want them to create graphs unto you”. This means that a good graph is a simple graph, in the Einsteinian sense that a graph should be made as simple as possible, but not simpler. A good graph communicates the main message effectively, without fuss and distraction. In addition, a good graph balances its graphical and textual elements – large symbols demand an increase in line width, and these together require an increase in font size.

The graphing chaos is exacerbated by the default settings in R (and the graphical packages that it provides, such as ggplot2), which are decidedly suboptimal. For instance, the font size is often too small, and the graphical elements are not sufficiently prominent. As a result, creating a good graph in R requires a lot of tinkering, not unlike the process of editing the first draft of a novice writer.

Fortunately, many plots share the same underlying structure, and the tinkering that has led to a clean graph of time series A will generally provide useful starting values for a clean graph of time series B. To exploit the overlap in structure, however, the user needs to remember the settings that were used for the first graph. Usually, this means that the user has to recall the location of the relevant R code. Sometimes the search for this initial code can take longer than the tinkering that was required to produce a clean graph in the first place.

In order to reduce the time needed to find relevant R code, we have constructed a compendium of clean graphs in R. This compendium, available at http://shinyapps.org/apps/RGraphCompendium/index.html, can also be used for teaching or as inspiration for improving one’s own graphs. In addition, the compendium provides a selective overview of the kind of graphs that researchers often use; the graphs cover a range of statistical scenarios and feature contributions of different data analysts. We do not wish to presume the graphs in the compendium are in any way perfect; some are better than others, and overall much remains to be improved. The compendium is undergoing continual refinement. Nevertheless, we hope the graphs are useful in their current state.

As an example of what the compendium has to offer, consider the graph below. This graph shows the proportion of the popular vote as a function of the relative height of the US president against his most successful opponent. Note the large circles for the data, the thick line for the linear relation, and the large font size for the axis labels. Also, note that the line does not touch the y-axis (a subtlety that requires deviating from the default). As in the compendium, the R code that created the graph is displayed after clicking the box “Show R-code”.

Show R-Code

# Presidential data up to and including 2008; data from Stulp et al. 2013

# rm(list=ls())

# height of president divided by height of most successful opponent:

height.ratio <- c(0.924324324, 1.081871345, 1, 0.971098266, 1.029761905,

0.935135135, 0.994252874, 0.908163265, 1.045714286, 1.18404908,

1.115606936, 0.971910112, 0.97752809, 0.978609626, 1,

0.933333333, 1.071428571, 0.944444444, 0.944444444, 1.017142857,

1.011111111, 1.011235955, 1.011235955, 1.089285714, 0.988888889,

1.011111111, 1.032967033, 1.044444444, 1, 1.086705202,

1.011560694, 1.005617978, 1.005617978, 1.005494505, 1.072222222,

1.011111111, 0.983783784, 0.967213115, 1.04519774, 1.027777778,

1.086705202, 1, 1.005347594, 0.983783784, 0.943005181, 1.057142857)

# proportion popular vote for president vs most successful opponent

# NB can be lower than .5 because popolar vote does not decide election

pop.vote <- c(0.427780852, 0.56148981, 0.597141922, 0.581254292, 0.530344067,

0.507425996, 0.526679292, 0.536690951, 0.577825976, 0.573225387,

0.550410082, 0.559380032, 0.484823958, 0.500466176, 0.502934212,

0.49569636, 0.516904414, 0.522050547, 0.531494442, 0.60014892,

0.545079801, 0.604274986, 0.51635906, 0.63850958, 0.652184407,

0.587920412, 0.5914898, 0.624614752, 0.550040193, 0.537771958,

0.523673642, 0.554517134, 0.577511576, 0.500856251, 0.613444534,

0.504063153, 0.617883695, 0.51049949, 0.553073235, 0.59166415,

0.538982024, 0.53455133, 0.547304058, 0.497350649, 0.512424242,

0.536914796)

#cor.test(height.ratio,pop.vote)

require(plotrix) # package plotrix is needed for function "ablineclip""

# if the following line and the line containing "dev.off()" are executed, the plot will be saved as a png file in the current working directory

# png("Presidental.png", width = 18, height = 18, units = "cm", res = 800, pointsize = 10)

op <- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5 , font.lab = 2, cex.axis = 1.3, bty = "n", las=1)

plot(height.ratio, pop.vote, col="black", pch=21, bg = "grey", cex = 2,

xlim=c(.90,1.20), ylim=c(.40,.70), ylab="", xlab="", axes=F)

axis(1)

axis(2)

reg1 <- lm(pop.vote~height.ratio)

ablineclip(reg1, lwd=2,x1 = .9, x2 = 1.2)

par(las=0)

mtext("Presidential Height Ratio", side=1, line=2.5, cex=1.5)

mtext("Relative Support for President", side=2, line=3.7, cex=1.5)

text(1.15, .65, "r = .39", cex=1.5)

# dev.off()

# For comparison, consider the default plot:

#par(op) # reset to default "par" settings

#plot(height.ratio, pop.vote) #yuk!

# rm(list=ls())

# height of president divided by height of most successful opponent:

height.ratio <- c(0.924324324, 1.081871345, 1, 0.971098266, 1.029761905,

0.935135135, 0.994252874, 0.908163265, 1.045714286, 1.18404908,

1.115606936, 0.971910112, 0.97752809, 0.978609626, 1,

0.933333333, 1.071428571, 0.944444444, 0.944444444, 1.017142857,

1.011111111, 1.011235955, 1.011235955, 1.089285714, 0.988888889,

1.011111111, 1.032967033, 1.044444444, 1, 1.086705202,

1.011560694, 1.005617978, 1.005617978, 1.005494505, 1.072222222,

1.011111111, 0.983783784, 0.967213115, 1.04519774, 1.027777778,

1.086705202, 1, 1.005347594, 0.983783784, 0.943005181, 1.057142857)

# proportion popular vote for president vs most successful opponent

# NB can be lower than .5 because popolar vote does not decide election

pop.vote <- c(0.427780852, 0.56148981, 0.597141922, 0.581254292, 0.530344067,

0.507425996, 0.526679292, 0.536690951, 0.577825976, 0.573225387,

0.550410082, 0.559380032, 0.484823958, 0.500466176, 0.502934212,

0.49569636, 0.516904414, 0.522050547, 0.531494442, 0.60014892,

0.545079801, 0.604274986, 0.51635906, 0.63850958, 0.652184407,

0.587920412, 0.5914898, 0.624614752, 0.550040193, 0.537771958,

0.523673642, 0.554517134, 0.577511576, 0.500856251, 0.613444534,

0.504063153, 0.617883695, 0.51049949, 0.553073235, 0.59166415,

0.538982024, 0.53455133, 0.547304058, 0.497350649, 0.512424242,

0.536914796)

#cor.test(height.ratio,pop.vote)

require(plotrix) # package plotrix is needed for function "ablineclip""

# if the following line and the line containing "dev.off()" are executed, the plot will be saved as a png file in the current working directory

# png("Presidental.png", width = 18, height = 18, units = "cm", res = 800, pointsize = 10)

op <- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5 , font.lab = 2, cex.axis = 1.3, bty = "n", las=1)

plot(height.ratio, pop.vote, col="black", pch=21, bg = "grey", cex = 2,

xlim=c(.90,1.20), ylim=c(.40,.70), ylab="", xlab="", axes=F)

axis(1)

axis(2)

reg1 <- lm(pop.vote~height.ratio)

ablineclip(reg1, lwd=2,x1 = .9, x2 = 1.2)

par(las=0)

mtext("Presidential Height Ratio", side=1, line=2.5, cex=1.5)

mtext("Relative Support for President", side=2, line=3.7, cex=1.5)

text(1.15, .65, "r = .39", cex=1.5)

# dev.off()

# For comparison, consider the default plot:

#par(op) # reset to default "par" settings

#plot(height.ratio, pop.vote) #yuk!

A more complicated example takes the same data, but uses it to plot the development of the Bayes factor, assessing the evidence for the hypothesis that taller presidential candidates attract more votes. This plot was created based in part on code from Ruud Wetzels and Benjamin Scheibehenne. Note the annotations on the right side of the plot, and the subtle horizontal lines that indicate Jeffreys’ criteria on the evidence. It took some time to figure out how to display the word “Evidence” in its current direction.

Show R-Code

# rm(list=ls())

# height of president divided by height of most successful opponent:

height.ratio <- c(0.924324324, 1.081871345, 1, 0.971098266, 1.029761905, 0.935135135, 0.994252874, 0.908163265, 1.045714286, 1.18404908, 1.115606936, 0.971910112, 0.97752809, 0.978609626, 1, 0.933333333, 1.071428571, 0.944444444, 0.944444444, 1.017142857, 1.011111111, 1.011235955, 1.011235955, 1.089285714, 0.988888889, 1.011111111, 1.032967033, 1.044444444, 1, 1.086705202, 1.011560694, 1.005617978, 1.005617978, 1.005494505, 1.072222222, 1.011111111, 0.983783784, 0.967213115, 1.04519774, 1.027777778, 1.086705202, 1, 1.005347594, 0.983783784, 0.943005181, 1.057142857)

# proportion popular vote for president vs most successful opponent

pop.vote <- c(0.427780852, 0.56148981, 0.597141922, 0.581254292, 0.530344067, 0.507425996, 0.526679292, 0.536690951, 0.577825976, 0.573225387, 0.550410082, 0.559380032, 0.484823958, 0.500466176, 0.502934212, 0.49569636, 0.516904414, 0.522050547, 0.531494442, 0.60014892, 0.545079801, 0.604274986, 0.51635906, 0.63850958, 0.652184407, 0.587920412, 0.5914898, 0.624614752, 0.550040193, 0.537771958, 0.523673642, 0.554517134, 0.577511576, 0.500856251, 0.613444534, 0.504063153, 0.617883695, 0.51049949, 0.553073235, 0.59166415, 0.538982024, 0.53455133, 0.547304058, 0.497350649, 0.512424242, 0.536914796)

## now calculate BF sequentially; two-sided test

library("hypergeo")

BF10.HG.exact = function(n, r)

{

#Jeffreys' test for whether a correlation is zero or not

#Jeffreys (1961), pp. 289-292

#Note that if the means are subtracted, n needs to be replaced by n-1

hypgeo = hypergeo((.25+n/2), (-.25+n/2), (3/2+n/2), r^2)

BF10 = ( sqrt(pi) * gamma(n/2+1) * (hypgeo) ) / ( 2 * gamma(3/2+n/2) )

return(as.numeric(BF10))

}

BF10 <- array()

BF10[1]<-1

BF10[2]<-1

for (i in 3:length(height.ratio))

{

BF10[i] <- BF10.HG.exact(n=i-1, r=cor(height.ratio[1:i],pop.vote[1:i]))

}

# We wish to plot this Bayes factor sequentially, as it unfolds as more elections become available:

#============ Plot log Bayes factors ===========================

par(cex.main = 1.3, mar = c(4.5, 6, 4, 7)+.1, mgp = c(3, 1, 0), #bottom, left, top, right

cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, las=1)

xhigh <- 60

plot(log(BF10), xlim=c(1,xhigh), ylim=c(-1*log(200),log(200)), xlab="", ylab="", cex.lab=1.3,cex.axis=1.3, las =1, yaxt="n", bty = "n", type="p", pch=21, bg="grey")

labelsUpper=log(c(100,30,10,3,1))

labelsLower=-1*labelsUpper

criticalP=c(labelsLower,0,labelsUpper)

for (idx in 1:length(criticalP))

{

abline(h=criticalP[idx],col='darkgrey',lwd=1,lty=2)

}

abline(h=0)

axis(side=4, at=criticalP,tick=T,las=2,cex.axis=1, labels=F)

axis(side=4, at=labelsUpper+.602, tick=F, cex.axis=1, labels=c("Extreme","Very strong", "Strong","Moderate", "Anecdotal"))

axis(side=4, at=labelsLower-.602,tick=F, cex.axis=1, labels=c("Extreme","Very strong", "Strong","Moderate", "Anecdotal"))

axis(side=2, at=c(criticalP),tick=T,las=2,cex.axis=1,

labels=c("1/100","1/30","1/10","1/3","1","", "100","30","10","3",""))

mtext(expression(BF[1][0]), side=2, line=2.5, las=0, cex=1.3)

grid::grid.text("Evidence", 0.97, 0.5, rot = 270, gp=grid::gpar(cex=1.3))

mtext("No. of Elections", side=1, line=2.5, las=1, cex=1.3)

arrows(20, -log(10), 20, -log(100), length=.25, angle=30, code=2, lwd=2)

arrows(20, log(10), 20, log(100), length=.25, angle=30, code=2, lwd=2)

text(25, -log(70), "Evidence for H0", pos=4, cex=1.3)

text(25, log(70), "Evidence for H1", pos=4, cex=1.3)

# height of president divided by height of most successful opponent:

height.ratio <- c(0.924324324, 1.081871345, 1, 0.971098266, 1.029761905, 0.935135135, 0.994252874, 0.908163265, 1.045714286, 1.18404908, 1.115606936, 0.971910112, 0.97752809, 0.978609626, 1, 0.933333333, 1.071428571, 0.944444444, 0.944444444, 1.017142857, 1.011111111, 1.011235955, 1.011235955, 1.089285714, 0.988888889, 1.011111111, 1.032967033, 1.044444444, 1, 1.086705202, 1.011560694, 1.005617978, 1.005617978, 1.005494505, 1.072222222, 1.011111111, 0.983783784, 0.967213115, 1.04519774, 1.027777778, 1.086705202, 1, 1.005347594, 0.983783784, 0.943005181, 1.057142857)

# proportion popular vote for president vs most successful opponent

pop.vote <- c(0.427780852, 0.56148981, 0.597141922, 0.581254292, 0.530344067, 0.507425996, 0.526679292, 0.536690951, 0.577825976, 0.573225387, 0.550410082, 0.559380032, 0.484823958, 0.500466176, 0.502934212, 0.49569636, 0.516904414, 0.522050547, 0.531494442, 0.60014892, 0.545079801, 0.604274986, 0.51635906, 0.63850958, 0.652184407, 0.587920412, 0.5914898, 0.624614752, 0.550040193, 0.537771958, 0.523673642, 0.554517134, 0.577511576, 0.500856251, 0.613444534, 0.504063153, 0.617883695, 0.51049949, 0.553073235, 0.59166415, 0.538982024, 0.53455133, 0.547304058, 0.497350649, 0.512424242, 0.536914796)

## now calculate BF sequentially; two-sided test

library("hypergeo")

BF10.HG.exact = function(n, r)

{

#Jeffreys' test for whether a correlation is zero or not

#Jeffreys (1961), pp. 289-292

#Note that if the means are subtracted, n needs to be replaced by n-1

hypgeo = hypergeo((.25+n/2), (-.25+n/2), (3/2+n/2), r^2)

BF10 = ( sqrt(pi) * gamma(n/2+1) * (hypgeo) ) / ( 2 * gamma(3/2+n/2) )

return(as.numeric(BF10))

}

BF10 <- array()

BF10[1]<-1

BF10[2]<-1

for (i in 3:length(height.ratio))

{

BF10[i] <- BF10.HG.exact(n=i-1, r=cor(height.ratio[1:i],pop.vote[1:i]))

}

# We wish to plot this Bayes factor sequentially, as it unfolds as more elections become available:

#============ Plot log Bayes factors ===========================

par(cex.main = 1.3, mar = c(4.5, 6, 4, 7)+.1, mgp = c(3, 1, 0), #bottom, left, top, right

cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, las=1)

xhigh <- 60

plot(log(BF10), xlim=c(1,xhigh), ylim=c(-1*log(200),log(200)), xlab="", ylab="", cex.lab=1.3,cex.axis=1.3, las =1, yaxt="n", bty = "n", type="p", pch=21, bg="grey")

labelsUpper=log(c(100,30,10,3,1))

labelsLower=-1*labelsUpper

criticalP=c(labelsLower,0,labelsUpper)

for (idx in 1:length(criticalP))

{

abline(h=criticalP[idx],col='darkgrey',lwd=1,lty=2)

}

abline(h=0)

axis(side=4, at=criticalP,tick=T,las=2,cex.axis=1, labels=F)

axis(side=4, at=labelsUpper+.602, tick=F, cex.axis=1, labels=c("Extreme","Very strong", "Strong","Moderate", "Anecdotal"))

axis(side=4, at=labelsLower-.602,tick=F, cex.axis=1, labels=c("Extreme","Very strong", "Strong","Moderate", "Anecdotal"))

axis(side=2, at=c(criticalP),tick=T,las=2,cex.axis=1,

labels=c("1/100","1/30","1/10","1/3","1","", "100","30","10","3",""))

mtext(expression(BF[1][0]), side=2, line=2.5, las=0, cex=1.3)

grid::grid.text("Evidence", 0.97, 0.5, rot = 270, gp=grid::gpar(cex=1.3))

mtext("No. of Elections", side=1, line=2.5, las=1, cex=1.3)

arrows(20, -log(10), 20, -log(100), length=.25, angle=30, code=2, lwd=2)

arrows(20, log(10), 20, log(100), length=.25, angle=30, code=2, lwd=2)

text(25, -log(70), "Evidence for H0", pos=4, cex=1.3)

text(25, log(70), "Evidence for H1", pos=4, cex=1.3)

A final example is borrowed from the graphs in **JASP** (http://jasp-stats.org), a free and open-source statistical software program with a GUI not unlike that of SPSS. In contrast to SPSS, JASP also includes Bayesian hypthesis tests, the results of which are summarized in graphs such as the one below.

Show R-Code

.plotPosterior.ttest <- function(x= NULL, y= NULL, paired= FALSE, oneSided= FALSE, iterations= 10000, rscale= "medium", lwd= 2, cexPoints= 1.5, cexAxis= 1.2, cexYlab= 1.5, cexXlab= 1.5, cexTextBF= 1.4, cexCI= 1.1, cexLegend= 1.4, lwdAxis= 1.2){

library(BayesFactor)

if(rscale == "medium"){

r <- sqrt(2) / 2

}

if(rscale == "wide"){

r <- 1

}

if(rscale == "ultrawide"){

r <- sqrt(2)

}

if(mode(rscale) == "numeric"){

r <- rscale

}

if(oneSided == FALSE){

nullInterval <- NULL

}

if(oneSided == "right"){

nullInterval <- c(0, Inf)

}

if(oneSided == "left"){

nullInterval <- c(-Inf, 0)

}

# sample from delta posterior

samples <- BayesFactor::ttestBF(x=x, y=y, paired=paired, nullInterval= nullInterval, posterior = TRUE, iterations = iterations, rscale= r)

delta <- samples[,"delta"]

# fit denisty estimator

fit.posterior <- logspline::logspline(delta)

# density function posterior

dposterior <- function(x, oneSided= oneSided, delta= delta){

if(oneSided == FALSE){

k <- 1

return(k*logspline::dlogspline(x, fit.posterior))

}

if(oneSided == "right"){

k <- 1 / (length(delta[delta >= 0]) / length(delta))

return(ifelse(x < 0, 0, k*logspline::dlogspline(x, fit.posterior)))

}

if(oneSided == "left"){

k <- 1 / (length(delta[delta <= 0]) / length(delta))

return(ifelse(x > 0, 0, k*logspline::dlogspline(x, fit.posterior)))

}

}

# pdf cauchy prior

dprior <- function(delta,r, oneSided= oneSided){

if(oneSided == "right"){

y <- ifelse(delta < 0, 0, 2/(pi*r*(1+(delta/r)^2)))

return(y)

}

if(oneSided == "left"){

y <- ifelse(delta > 0, 0, 2/(pi*r*(1+(delta/r)^2)))

return(y)

} else{

return(1/(pi*r*(1+(delta/r)^2)))

}

}

# set limits plot

xlim <- vector("numeric", 2)

if(oneSided == FALSE){

xlim[1] <- min(-2, quantile(delta, probs = 0.01)[[1]])

xlim[2] <- max(2, quantile(delta, probs = 0.99)[[1]])

}

if(oneSided == "right"){

xlim[1] <- min(-2, quantile(delta[delta >= 0], probs = 0.01)[[1]])

xlim[2] <- max(2, quantile(delta[delta >= 0], probs = 0.99)[[1]])

}

if(oneSided == "left"){

xlim[1] <- min(-2, quantile(delta[delta <= 0], probs = 0.01)[[1]])

xlim[2] <- max(2, quantile(delta[delta <= 0], probs = 0.99)[[1]])

}

ylim <- vector("numeric", 2)

ylim[1] <- 0

ylim[2] <- max(dprior(0,r, oneSided= oneSided), 1.28*max(dposterior(x= delta, oneSided= oneSided, delta=delta)))

# calculate position of "nice" tick marks and create labels

xticks <- pretty(xlim)

yticks <- pretty(ylim)

xlabels <- formatC(pretty(xlim), 1, format= "f")

ylabels <- formatC(pretty(ylim), 1, format= "f")

# 95% credible interval:

if(oneSided == FALSE){

CIlow <- quantile(delta, probs = 0.025)[[1]]

CIhigh <- quantile(delta, probs = 0.975)[[1]]

}

if(oneSided == "right"){

CIlow <- quantile(delta[delta >= 0], probs = 0.025)[[1]]

CIhigh <- quantile(delta[delta >= 0], probs = 0.975)[[1]]

}

if(oneSided == "left"){

CIlow <- quantile(delta[delta <= 0], probs = 0.025)[[1]]

CIhigh <- quantile(delta[delta <= 0], probs = 0.975)[[1]]

}

par(mar= c(5, 5, 7, 4) + 0.1, las=1)

xlim <- c(min(CIlow,range(xticks)[1]), max(range(xticks)[2], CIhigh))

plot(1,1, xlim= xlim, ylim= range(yticks), ylab= "", xlab="", type= "n", axes= FALSE)

lines(seq(min(xticks), max(xticks),length.out = 1000),dposterior(x=seq(min(xticks), max(xticks),length.out = 1000), oneSided = oneSided, delta=delta), lwd= lwd, xlim= xlim, ylim= range(yticks), ylab= "", xlab= "")

lines(seq(min(xticks), max(xticks),length.out = 1000), dprior(seq(min(xticks), max(xticks),length.out = 1000), r=r, oneSided= oneSided), lwd= lwd, lty=3)

axis(1, at= xticks, labels = xlabels, cex.axis= cexAxis, lwd= lwdAxis)

axis(2, at= yticks, labels= ylabels, , cex.axis= cexAxis, lwd= lwdAxis)

mtext(text = "Density", side = 2, las=0, cex = cexYlab, line= 3)

mtext(expression(paste("Effect size", ~delta)), side = 1, cex = cexXlab, line= 2.5)

points(0, dprior(0,r, oneSided= oneSided), col="black", pch=21, bg = "grey", cex= cexPoints)

points(0, dposterior(0, oneSided = oneSided, delta=delta), col="black", pch=21, bg = "grey", cex= cexPoints)

# 95% credible interval

dmax <- optimize(function(x)dposterior(x,oneSided= oneSided, delta=delta), interval= range(xticks), maximum = TRUE)$objective # get maximum density

yCI <- grconvertY(dmax, "user", "ndc") + 0.08

yCIt <- grconvertY(dmax, "user", "ndc") + 0.04

y95 <- grconvertY(dmax, "user", "ndc") + 0.1

yCI <- grconvertY(yCI, "ndc", "user")

yCIt <- grconvertY(yCIt, "ndc", "user")

y95 <- grconvertY(y95, "ndc", "user")

arrows(CIlow, yCI , CIhigh, yCI, angle = 90, code = 3, length= 0.1, lwd= lwd)

text(mean(c(CIlow, CIhigh)), y95,"95%", cex= cexCI)

text(CIlow, yCIt, bquote(.(formatC(CIlow,2, format="f"))), cex= cexCI)

text(CIhigh, yCIt, bquote(.(formatC(CIhigh,2, format= "f"))), cex= cexCI)

# enable plotting in margin

par(xpd=TRUE)

# display BF10 value

BF <- BayesFactor::ttestBF(x=x, y=y, paired=paired, nullInterval= nullInterval, posterior = FALSE, rscale= r)

BF10 <- BayesFactor::extractBF(BF, logbf = FALSE, onlybf = F)[1, "bf"]

BF01 <- 1 / BF10

xx <- grconvertX(0.3, "ndc", "user")

yy <- grconvertY(0.822, "ndc", "user")

yy2 <- grconvertY(0.878, "ndc", "user")

if(BF10 >= 1000000 | BF01 >= 1000000){

BF10t <- format(BF10, digits= 3, scientific = TRUE)

BF01t <- format(BF01, digits= 3, scientific = TRUE)

}

if(BF10 < 1000000 & BF01 < 1000000){

BF10t <- formatC(BF10,2, format = "f")

BF01t <- formatC(BF01,2, format = "f")

}

if(oneSided == FALSE){

text(xx, yy2, bquote(BF[10]==.(BF10t)), cex= cexTextBF)

text(xx, yy, bquote(BF[0][1]==.(BF01t)), cex= cexTextBF)

}

if(oneSided == "right"){

text(xx, yy2, bquote(BF["+"][0]==.(BF10t)), cex= cexTextBF)

text(xx, yy, bquote(BF[0]["+"]==.(BF01t)), cex= cexTextBF)

}

if(oneSided == "left"){

text(xx, yy2, bquote(BF["-"][0]==.(BF10t)), cex= cexTextBF)

text(xx, yy, bquote(BF[0]["-"]==.(BF01t)), cex= cexTextBF)

}

# probability wheel

if(max(nchar(BF10t), nchar(BF01t)) <= 4){

xx <- grconvertX(0.44, "ndc", "user")

}

# probability wheel

if(max(nchar(BF10t), nchar(BF01t)) == 5){

xx <- grconvertX(0.44 + 0.001* 5, "ndc", "user")

}

# probability wheel

if(max(nchar(BF10t), nchar(BF01t)) == 6){

xx <- grconvertX(0.44 + 0.001* 6, "ndc", "user")

}

if(max(nchar(BF10t), nchar(BF01t)) == 7){

xx <- grconvertX(0.44 + 0.002* max(nchar(BF10t), nchar(BF01t)), "ndc", "user")

}

if(max(nchar(BF10t), nchar(BF01t)) == 8){

xx <- grconvertX(0.44 + 0.003* max(nchar(BF10t), nchar(BF01t)), "ndc", "user")

}

if(max(nchar(BF10t), nchar(BF01t)) > 8){

xx <- grconvertX(0.44 + 0.004* max(nchar(BF10t), nchar(BF01t)), "ndc", "user")

}

yy <- grconvertY(0.85, "ndc", "user")

# make sure that colored area is centered

radius <- 0.06*diff(range(xticks))

A <- radius^2*pi

alpha <- 2 / (BF01 + 1) * A / radius^2

startpos <- pi/2 - alpha/2

# draw probability wheel

plotrix::floating.pie(xx, yy,c(BF10, 1),radius= radius, col=c("darkred", "white"), lwd=2,startpos = startpos)

yy <- grconvertY(0.927, "ndc", "user")

yy2 <- grconvertY(0.77, "ndc", "user")

if(oneSided == FALSE){

text(xx, yy, "data|H1", cex= cexCI)

text(xx, yy2, "data|H0", cex= cexCI)

}

if(oneSided == "right"){

text(xx, yy, "data|H+", cex= cexCI)

text(xx, yy2, "data|H0", cex= cexCI)

}

if(oneSided == "left"){

text(xx, yy, "data|H-", cex= cexCI)

text(xx, yy2, "data|H0", cex= cexCI)

}

# add legend

xx <- grconvertX(0.57, "ndc", "user")

yy <- grconvertY(0.92, "ndc", "user")

legend(xx, yy, legend = c("Posterior", "Prior"), lty=c(1,3), bty= "n", lwd = c(lwd,lwd), cex= cexLegend)

}

set.seed(1)

.plotPosterior.ttest(x= rnorm(30,0.15), rscale=1)

library(BayesFactor)

if(rscale == "medium"){

r <- sqrt(2) / 2

}

if(rscale == "wide"){

r <- 1

}

if(rscale == "ultrawide"){

r <- sqrt(2)

}

if(mode(rscale) == "numeric"){

r <- rscale

}

if(oneSided == FALSE){

nullInterval <- NULL

}

if(oneSided == "right"){

nullInterval <- c(0, Inf)

}

if(oneSided == "left"){

nullInterval <- c(-Inf, 0)

}

# sample from delta posterior

samples <- BayesFactor::ttestBF(x=x, y=y, paired=paired, nullInterval= nullInterval, posterior = TRUE, iterations = iterations, rscale= r)

delta <- samples[,"delta"]

# fit denisty estimator

fit.posterior <- logspline::logspline(delta)

# density function posterior

dposterior <- function(x, oneSided= oneSided, delta= delta){

if(oneSided == FALSE){

k <- 1

return(k*logspline::dlogspline(x, fit.posterior))

}

if(oneSided == "right"){

k <- 1 / (length(delta[delta >= 0]) / length(delta))

return(ifelse(x < 0, 0, k*logspline::dlogspline(x, fit.posterior)))

}

if(oneSided == "left"){

k <- 1 / (length(delta[delta <= 0]) / length(delta))

return(ifelse(x > 0, 0, k*logspline::dlogspline(x, fit.posterior)))

}

}

# pdf cauchy prior

dprior <- function(delta,r, oneSided= oneSided){

if(oneSided == "right"){

y <- ifelse(delta < 0, 0, 2/(pi*r*(1+(delta/r)^2)))

return(y)

}

if(oneSided == "left"){

y <- ifelse(delta > 0, 0, 2/(pi*r*(1+(delta/r)^2)))

return(y)

} else{

return(1/(pi*r*(1+(delta/r)^2)))

}

}

# set limits plot

xlim <- vector("numeric", 2)

if(oneSided == FALSE){

xlim[1] <- min(-2, quantile(delta, probs = 0.01)[[1]])

xlim[2] <- max(2, quantile(delta, probs = 0.99)[[1]])

}

if(oneSided == "right"){

xlim[1] <- min(-2, quantile(delta[delta >= 0], probs = 0.01)[[1]])

xlim[2] <- max(2, quantile(delta[delta >= 0], probs = 0.99)[[1]])

}

if(oneSided == "left"){

xlim[1] <- min(-2, quantile(delta[delta <= 0], probs = 0.01)[[1]])

xlim[2] <- max(2, quantile(delta[delta <= 0], probs = 0.99)[[1]])

}

ylim <- vector("numeric", 2)

ylim[1] <- 0

ylim[2] <- max(dprior(0,r, oneSided= oneSided), 1.28*max(dposterior(x= delta, oneSided= oneSided, delta=delta)))

# calculate position of "nice" tick marks and create labels

xticks <- pretty(xlim)

yticks <- pretty(ylim)

xlabels <- formatC(pretty(xlim), 1, format= "f")

ylabels <- formatC(pretty(ylim), 1, format= "f")

# 95% credible interval:

if(oneSided == FALSE){

CIlow <- quantile(delta, probs = 0.025)[[1]]

CIhigh <- quantile(delta, probs = 0.975)[[1]]

}

if(oneSided == "right"){

CIlow <- quantile(delta[delta >= 0], probs = 0.025)[[1]]

CIhigh <- quantile(delta[delta >= 0], probs = 0.975)[[1]]

}

if(oneSided == "left"){

CIlow <- quantile(delta[delta <= 0], probs = 0.025)[[1]]

CIhigh <- quantile(delta[delta <= 0], probs = 0.975)[[1]]

}

par(mar= c(5, 5, 7, 4) + 0.1, las=1)

xlim <- c(min(CIlow,range(xticks)[1]), max(range(xticks)[2], CIhigh))

plot(1,1, xlim= xlim, ylim= range(yticks), ylab= "", xlab="", type= "n", axes= FALSE)

lines(seq(min(xticks), max(xticks),length.out = 1000),dposterior(x=seq(min(xticks), max(xticks),length.out = 1000), oneSided = oneSided, delta=delta), lwd= lwd, xlim= xlim, ylim= range(yticks), ylab= "", xlab= "")

lines(seq(min(xticks), max(xticks),length.out = 1000), dprior(seq(min(xticks), max(xticks),length.out = 1000), r=r, oneSided= oneSided), lwd= lwd, lty=3)

axis(1, at= xticks, labels = xlabels, cex.axis= cexAxis, lwd= lwdAxis)

axis(2, at= yticks, labels= ylabels, , cex.axis= cexAxis, lwd= lwdAxis)

mtext(text = "Density", side = 2, las=0, cex = cexYlab, line= 3)

mtext(expression(paste("Effect size", ~delta)), side = 1, cex = cexXlab, line= 2.5)

points(0, dprior(0,r, oneSided= oneSided), col="black", pch=21, bg = "grey", cex= cexPoints)

points(0, dposterior(0, oneSided = oneSided, delta=delta), col="black", pch=21, bg = "grey", cex= cexPoints)

# 95% credible interval

dmax <- optimize(function(x)dposterior(x,oneSided= oneSided, delta=delta), interval= range(xticks), maximum = TRUE)$objective # get maximum density

yCI <- grconvertY(dmax, "user", "ndc") + 0.08

yCIt <- grconvertY(dmax, "user", "ndc") + 0.04

y95 <- grconvertY(dmax, "user", "ndc") + 0.1

yCI <- grconvertY(yCI, "ndc", "user")

yCIt <- grconvertY(yCIt, "ndc", "user")

y95 <- grconvertY(y95, "ndc", "user")

arrows(CIlow, yCI , CIhigh, yCI, angle = 90, code = 3, length= 0.1, lwd= lwd)

text(mean(c(CIlow, CIhigh)), y95,"95%", cex= cexCI)

text(CIlow, yCIt, bquote(.(formatC(CIlow,2, format="f"))), cex= cexCI)

text(CIhigh, yCIt, bquote(.(formatC(CIhigh,2, format= "f"))), cex= cexCI)

# enable plotting in margin

par(xpd=TRUE)

# display BF10 value

BF <- BayesFactor::ttestBF(x=x, y=y, paired=paired, nullInterval= nullInterval, posterior = FALSE, rscale= r)

BF10 <- BayesFactor::extractBF(BF, logbf = FALSE, onlybf = F)[1, "bf"]

BF01 <- 1 / BF10

xx <- grconvertX(0.3, "ndc", "user")

yy <- grconvertY(0.822, "ndc", "user")

yy2 <- grconvertY(0.878, "ndc", "user")

if(BF10 >= 1000000 | BF01 >= 1000000){

BF10t <- format(BF10, digits= 3, scientific = TRUE)

BF01t <- format(BF01, digits= 3, scientific = TRUE)

}

if(BF10 < 1000000 & BF01 < 1000000){

BF10t <- formatC(BF10,2, format = "f")

BF01t <- formatC(BF01,2, format = "f")

}

if(oneSided == FALSE){

text(xx, yy2, bquote(BF[10]==.(BF10t)), cex= cexTextBF)

text(xx, yy, bquote(BF[0][1]==.(BF01t)), cex= cexTextBF)

}

if(oneSided == "right"){

text(xx, yy2, bquote(BF["+"][0]==.(BF10t)), cex= cexTextBF)

text(xx, yy, bquote(BF[0]["+"]==.(BF01t)), cex= cexTextBF)

}

if(oneSided == "left"){

text(xx, yy2, bquote(BF["-"][0]==.(BF10t)), cex= cexTextBF)

text(xx, yy, bquote(BF[0]["-"]==.(BF01t)), cex= cexTextBF)

}

# probability wheel

if(max(nchar(BF10t), nchar(BF01t)) <= 4){

xx <- grconvertX(0.44, "ndc", "user")

}

# probability wheel

if(max(nchar(BF10t), nchar(BF01t)) == 5){

xx <- grconvertX(0.44 + 0.001* 5, "ndc", "user")

}

# probability wheel

if(max(nchar(BF10t), nchar(BF01t)) == 6){

xx <- grconvertX(0.44 + 0.001* 6, "ndc", "user")

}

if(max(nchar(BF10t), nchar(BF01t)) == 7){

xx <- grconvertX(0.44 + 0.002* max(nchar(BF10t), nchar(BF01t)), "ndc", "user")

}

if(max(nchar(BF10t), nchar(BF01t)) == 8){

xx <- grconvertX(0.44 + 0.003* max(nchar(BF10t), nchar(BF01t)), "ndc", "user")

}

if(max(nchar(BF10t), nchar(BF01t)) > 8){

xx <- grconvertX(0.44 + 0.004* max(nchar(BF10t), nchar(BF01t)), "ndc", "user")

}

yy <- grconvertY(0.85, "ndc", "user")

# make sure that colored area is centered

radius <- 0.06*diff(range(xticks))

A <- radius^2*pi

alpha <- 2 / (BF01 + 1) * A / radius^2

startpos <- pi/2 - alpha/2

# draw probability wheel

plotrix::floating.pie(xx, yy,c(BF10, 1),radius= radius, col=c("darkred", "white"), lwd=2,startpos = startpos)

yy <- grconvertY(0.927, "ndc", "user")

yy2 <- grconvertY(0.77, "ndc", "user")

if(oneSided == FALSE){

text(xx, yy, "data|H1", cex= cexCI)

text(xx, yy2, "data|H0", cex= cexCI)

}

if(oneSided == "right"){

text(xx, yy, "data|H+", cex= cexCI)

text(xx, yy2, "data|H0", cex= cexCI)

}

if(oneSided == "left"){

text(xx, yy, "data|H-", cex= cexCI)

text(xx, yy2, "data|H0", cex= cexCI)

}

# add legend

xx <- grconvertX(0.57, "ndc", "user")

yy <- grconvertY(0.92, "ndc", "user")

legend(xx, yy, legend = c("Posterior", "Prior"), lty=c(1,3), bty= "n", lwd = c(lwd,lwd), cex= cexLegend)

}

set.seed(1)

.plotPosterior.ttest(x= rnorm(30,0.15), rscale=1)

The compendium contains many more examples. We hope some R users will find them convenient. Finally, if you create a clean graph in R that you believe is a candidate for inclusion in this compendium, please do not hesitate to write an email to EJ.Wagenmakers@gmail.com. Your contribution will be acknowledged explicitly, alongside the code you provided.

*Eric-Jan Wagenmakers and Quentin Gronau*

*University of Amsterdam, Department of Psychology.*

Key links:

[…] A Compendium of Clean Graphs in R Every data analyst knows that a good graph is worth a thousand words, and perhaps a hundred tables. But how should one create a good, clean graph? In R, this task is anything but easy. Many users find it almost impossible to resist the siren song of adding grid lines, including grey backgrounds, using elaborate color schemes, and applying default font sizes that makes the text much too small in relation to the graphical elements. As a result, many R graphs are an aesthetic disaster; they are difficult to parse and unfit for publication. […]

Hello,

astonishingly nice plots. Is there any possibility to download the shiny app?

Actually, this is not a shiny app, but a rmarkdown file which has been compiled to html. All R source code is shown when you click on the the “Show R-Code” box.

I can’t contribute since I’m using python and matplotlib. However, I think it is a bad idea to omit grid from the plot. For instance, I guess that the authors of the scatter plot would not report the estimated quantity that is of most interest to me – by how many percent does the pres. support increase when the height ratio increases by one percent. It’s more difficult to compute this without the grid.

Very nice. You’ve done a ton of work. It is much appreciated. I too am a proponent of R graphics. The control over the primitives, placement, etc. is very important.

I’ve been using different ways of presenting R methods to students and colleagues with Rmarkdown, such as step-by-step building of plots over named chunks to aid pedagogy. I’ve been using a custom.css file to help with fonts, colors, etc. Petr Keil’s “Simple Template for Scientific Manuscripts in Rmarkdown” (http://www.petrkeil.com/?p=2401) has been useful.

I’m curious, how tedious was it to insert the “Show Code” javascript into the RMD file? I’d love to incorporate that into my own work.

Thanks,

Edward C.

University of Minnesota, U.S.A.

Oops, sorry matus–not a reply to you.

Edward C.

> I’m curious, how tedious was it to insert the “Show Code” javascript into the RMD file? I’d love to incorporate that into my own work.

Well, it was quite tedious until I figured out how to do it …

But now it works! Here’s the relevant snippet:

https://gist.github.com/nicebread/bb0979db92738303cdb0

I nearly got there with “View Source,” but this is better.

Thanks for sharing it,

Edward C.

I got this site from my pal who informed me about this site and at the moment this time I

am browsing this website and reading very informative

articles here.

my blog – tiles

[…] in Political Science Survey Experiments: Comparing Questionnaires to Published Results – A Compendium of Clean Graphs in R – Bringing the Gold Standard into the Classroom: Replication in University Teaching – Awesome […]

Very nice! Many thanks for this! I notice most of the code is quite verbose though, and mainly uses base R plot commands. Could be nice to do something equivalent, but using lattice and/or ggplot2. Especially with the themes of ggplot2 and ggthemes one can reasonably easily get something that looks nice to the eye (not least to get rid of the funny weird looking grey background, using + theme_few() ).