8

Can ggplot2 be used to produce a so-called topoplot (often used in neuroscience)?

topoplot

Sample data:

   label          x          y     signal
1     R3 0.64924459 0.91228430  2.0261520
2     R4 0.78789621 0.78234410  1.7880972
3     R5 0.93169511 0.72980685  0.9170998
4     R6 0.48406513 0.82383895  3.1933129

Full sample data.

Rows represent individual electrodes. Columns x and y represent the projection into 2D space and the column signal is essentially the z-axis representing voltage measured at a given electrode.

stat_contour doesn't work, apparently due to unequal grid.

geom_density_2d only provides a density estimation of x and y.

geom_raster is one not fitted for this task or I must be using it incorrectly since it quickly runs out of memory.

Smoothing (like in the image on the right) and head contours (nose, ears) aren't necessary.

I want to avoid Matlab and transforming the data so that it fits this or that toolbox… Many thanks!

Update (26 January 2016)

The closest I've been able to get to my objective is via

library(colorRamps)
ggplot(channels, aes(x, y, z = signal)) + stat_summary_2d() + scale_fill_gradientn(colours=matlab.like(20))

which produces an image like this:

enter image description here

Update 2 (27 January 2016)

I've tried @alexforrence's approach with full data and this is the result:

@alexforrence's approach

It's a great start but there is a couple of issues:

  1. The last call (ggplot()) takes about 40 seconds on an Intel i7 4790K while Matlab toolboxes manage to generate these almost instantly; my ‘emergency solution’ above takes about a second.
  2. As you can see, the upper and lower border of the central part appear to be ‘sliced’ – I'm not sure what causes this but it could be the third issue.
  3. I'm getting these warnings:

    1: Removed 170235 rows containing non-finite values (stat_contour). 
    2: Removed 170235 rows containing non-finite values (stat_contour). 
    

Update 3 (27 January 2016)

Comparison between two plots produced with different interp(xo, yo) and stat_contour(binwidth) values:

comparison between different values

Ragged edges if one chooses low interp(xo, yo), in this case xo/yo = seq(0, 1, length = 100):

ragged edges

epo3
  • 2,991
  • 2
  • 33
  • 60
Harold Cavendish
  • 869
  • 10
  • 22

1 Answers1

8

Here's a potential start:

First, we'll attach some packages. I'm using akima to do linear interpolation, though it looks like EEGLAB uses some sort of spherical interpolation here? (the data was a little sparse to try it).

library(ggplot2)
library(akima)
library(reshape2)

Next, reading in the data:

dat <- read.table(text = "   label          x          y     signal
1     R3 0.64924459 0.91228430  2.0261520
2     R4 0.78789621 0.78234410  1.7880972
3     R5 0.93169511 0.72980685  0.9170998
4     R6 0.48406513 0.82383895  3.1933129")

We'll interpolate the data, and stick that in a data frame.

datmat <- interp(dat$x, dat$y, dat$signal, 
                 xo = seq(0, 1, length = 1000),
                 yo = seq(0, 1, length = 1000))
datmat2 <- melt(datmat$z)
names(datmat2) <- c('x', 'y', 'value')
datmat2[,1:2] <- datmat2[,1:2]/1000 # scale it back

I'm going to borrow from some previous answers. The circleFun below is from Draw a circle with ggplot2.

circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
  r = diameter / 2
  tt <- seq(0,2*pi,length.out = npoints)
  xx <- center[1] + r * cos(tt)
  yy <- center[2] + r * sin(tt)
  return(data.frame(x = xx, y = yy))
}

circledat <- circleFun(c(.5, .5), 1, npoints = 100) # center on [.5, .5]

# ignore anything outside the circle
datmat2$incircle <- (datmat2$x - .5)^2 + (datmat2$y - .5)^2 < .5^2 # mark
datmat2 <- datmat2[datmat2$incircle,]

And I really liked the look of the contour plot in R plot filled.contour() output in ggpplot2, so we'll borrow that one.

ggplot(datmat2, aes(x, y, z = value)) +
  geom_tile(aes(fill = value)) +
  stat_contour(aes(fill = ..level..), geom = 'polygon', binwidth = 0.01) +
  geom_contour(colour = 'white', alpha = 0.5) +
  scale_fill_distiller(palette = "Spectral", na.value = NA) + 
  geom_path(data = circledat, aes(x, y, z = NULL)) +
  # draw the nose (haven't drawn ears yet)
  geom_line(data = data.frame(x = c(0.45, 0.5, .55), y = c(1, 1.05, 1)), 
            aes(x, y, z = NULL)) +
  # add points for the electrodes
  geom_point(data = dat, aes(x, y, z = NULL, fill = NULL), 
             shape = 21, colour = 'black', fill = 'white', size = 2) +
  theme_bw()

enter image description here


With improvements mentioned in the comments (setting extrap = TRUE and linear = FALSE in the interp call to fill in gaps and do a spline smoothing, respectively, and removing NAs before plotting), we get:

enter image description here


mgcv can do spherical splines. This replaces akima (the chunk containing interp() isn't necessary).

library(mgcv)
spl1 <- gam(signal ~ s(x, y, bs = 'sos'), data = dat)
# fine grid, coarser is faster
datmat2 <- data.frame(expand.grid(x = seq(0, 1, 0.001), y = seq(0, 1, 0.001)))
resp <- predict(spl1, datmat2, type = "response")
datmat2$value <- resp

enter image description here

Community
  • 1
  • 1
alexforrence
  • 2,724
  • 2
  • 27
  • 30
  • 1
    This is great, especially as a ‘potential start’, many thanks! Unfortunately, there seem to be a couple of issues, mainly that a) the last call takes about 40 seconds on an Intel i7 4790K while Matlab toolboxes manage to generate these almost instantly (why is that?), b) I'm getting warnings about removed rows, and c) the plot appears to be ‘sliced’ around the central part (see my updated question in a couple of minutes, please). In case you're willing to play with it, I've also added full data for a given time point just below the original sample. – Harold Cavendish Jan 27 '16 at 13:34
  • Brilliant, thank you! Not only it looks great and does the job but your description is also very educational. It still takes plenty of time but I've noticed that if I reduce `interp() xo/yo = length` to 100, it looks _practically_ the same and only takes about 3 seconds. Raising `stat_countour() binwidth to 0.1` makes it almost instantaneous. The edges are a bit rough (the bigger the image, the more ‘pixelated’ edges) but perhaps there's a way of extending the fill a bit outside the circle and then cutting it off by masking the outside? See my update. – Harold Cavendish Jan 27 '16 at 16:42
  • Also if you have the chance to look at how EEGLAB does it, it will be much appreciated. – Harold Cavendish Jan 27 '16 at 16:44
  • One more thing, if I may: at certain other time points of my data, the plot appears to contain patterns such as prominent straight lines between colour levels (I mean the contrast between colours) which shouldn't be there (even without the tweaks to speed it up). Is this due to issues with interpolation? Can it be fixed? – Harold Cavendish Jan 27 '16 at 20:31
  • Try the new splines and see what you get. The rough edges can be a new (separate) question. I'm going to delete some comments to clean up too. – alexforrence Jan 28 '16 at 16:30
  • I've been trying to get this to work but I'm getting an error at `te(bs = 'sos')`: `Error in smooth.construct.sos.smooth.spec(object$margin[[i]], dat, knt) : Can only deal with a sphere`. Most other smooth terms are working fine but I'm not going to pretend I understand what the differences are. The number of those straight contrast lines has been reduced though. – Harold Cavendish Jan 30 '16 at 18:15
  • 1
    I might have accidentally copied a failed tensor spline attempt, give it a go with `s()` rather than `te()`. – alexforrence Jan 30 '16 at 19:43
  • It works, thank you! You've been very helpful, I really appreciate your patience with this. – Harold Cavendish Jan 30 '16 at 21:41