The end of the line for error bars in R

When plotting in R, I often use the segments function to add lines representing confidence intervals. This is a very simple way to plot lines connecting pairs of x,y coordinates.

Recently I discovered that by default, segments are styled with rounded line caps, which add to their length. This means, of course, that confidence intervals are slightly wider than intended.

R provides three styles of line ending – round, butt and square – which can be specified by the lend argument. The figure here shows the outcome of using each line ending, with vertical lines indicating actual end-points of segments. Both round and square line ends overshoot these points, while butt ends represent them correctly.

plot.new()
par(mar=c(1, 4, 1, 1))
plot.window(xlim=c(0, 1), ylim=c(0.5, 3.5))
axis(2, 1:3, c('round', 'butt', 'square'), las=1)
box(lwd=2)
segments(0.1, 1, 0.9, 1, lwd=20, lend='round')
segments(0.1, 2, 0.9, 2, lwd=20, lend='butt')
segments(0.1, 3, 0.9, 3, lwd=20, lend='square')
abline(v=c(0.1, 0.9))
Line end styles applied to segments plotted in R. Only 'butt' accurately represents end points.

Line end styles applied to segments plotted in R. Only ‘butt’ accurately represents end points.

The effect is slight, and is emphasized when line width is large. Regardless, it’s a good idea to routinely add lend='butt' (or lend=1) to your segments function calls.

A secondary benefit is that lines will appear crisper than when plotted with the default round caps.

Getting rasters into shape from R

Update: Windows people should use the modified version of the function provided by Francisco Rodriguez-Sanchez (mentioned in this comment), which is available as a Gist here.


Today I needed to convert a raster to a polygon shapefile for further processing and plotting in R. I like to keep my code together so I can easily keep track of what I’ve done, so it made sense to do the conversion in R as well. The fantastic raster package has a function that can do this (rasterToPolygons), but it seems to take a very, very long time and chew up all system resources (at least on my not particularly souped up X200 laptop).

We can cut down the time taken for conversion by calling a GDAL utility,  gdal_polygonize.py, directly from R using system2(). GDAL needs to be installed on your system, but you’ll probably want it installed anyway if you’re planning to talk to spatial data (in fact, you might find you actually already have it installed, as it is required for a bunch of other GIS software packages). Together with it’s cousin OGR, GDAL allows reading and writing of a wide range of raster and vector geospatial data formats. Once installed, you have at your disposal a bunch of GDAL Utilities that can be run from the terminal, including gdal_polygonize.py. For polygonize and a couple of others, you’ll also need Python installed (again, you may already have it installed, so check first!). Finally, it helps if the path to the gdal_polygonize.py exists in the global PATH variable. To check if this is the case, run the following in R:

Sys.which("gdal_polygonize.py")

If this returns an empty string, it looks like gdal_polygonize.py either doesn’t exist on your system, or the path to its containing directoy is missing from the PATH variable. In the latter case, you can either use the pypath argument accepted by the function below to specify the path, or modify your PATH variable. The PATH variable can be modified by following these instructions for Windows and Mac. I suspect Linux users should typically not run into this problem, as the GDAL executables are usually installed to paths that already exist in PATH (e.g. /usr/bin).

Anyhow… let’s compare rasterToPolygons with gdal_polygonize. The function included below (gdal_polygonizeR) borrows from code provided in a post by Lyndon Estes on R-sig-geo. Thanks mate!

Continue reading

R functions to filter rjags results

A while back I was running a bunch of JAGS models through R, using the rjags (written by Martyn Plummer) and R2jags (by Yu-Sung Su) packages. These packages provide a great interface to the JAGS software, which allows analysis of Bayesian models (written in the BUGS language) through Markov chain Monte Carlo simulation.

Running a JAGS model using these tools returns an rjags object, which when printed to the screen, summarises the posterior distribution of each monitored node, giving its mean and standard deviation, a range of quantiles, and its Gelman-Rubin convergence diagnostic statistic (Rhat), which indicates the ratio of variance within chains to that among chains. The summary is great, but when monitoring a large number of nodes, printing these to the screen can cause R to hang, and can exceed the screen buffer (not to mention making it painful to find the nodes you’re immediately interested in).

To help deal with this I wrote a couple of simple R functions:

  • jagsresults: return a matrix containing summary results for just the nodes you are interested in (using regular expression pattern-matching, if desired).
  • rhats: sort the output by the nodes’ Rhat values, making it easy to show the n least converged nodes.

Continue reading