How to plot a peristimulus time histogram (PSTH) in R with ggplot2

Question

Say that I have two conditions, 'a' and 'b'. A neuron fires on average 40 spikes / second (Hz) in condition 'a' and 80 spikes / second in condition 'b'. The response to condition 'a' is presented 20 times and condition 'b' is presented 10 times, with each presentation for 1000 ms.

AB <- rbind(
    ldply( 1:20, 
        function(trial) { 
          data.frame( 
              trial=trial, 
              cond=factor('a',c('a','b')), 
              spiketime = runif(40,0,1000))
        }
    ), ldply(21:30, 
        function(trial) {
          data.frame(
              trial=trial, 
              cond=factor('b',c('a','b')), 
              spiketime = runif(80,0,1000))
        }
  )
)

A simple histogram can be plotted with:

qplot(spiketime, data=AB, geom='line',stat='bin',y=..count.., 
      xlim=c(0,1000), colour=cond, binwidth=100,xlab='Milliseconds')

However, this does not average over the presentations, and hence, the values on the y axis are about the same. I would like to plot the spiking rate (spikes/second) along the y-axis, which would show that condition 'b' elicits about twice as many spikes per second. The spiking rate does not increase as number of presentations increase, it just becomes less noisy. Is there a way to do this without pre-processing the dataframe AB?

In other words, can I do something along the lines of:

qplot(spiketime, data=AB, geom='line',stat='bin',
      y=..count../num_presentations*1000/binwidth, ylab='Spikes per second',
      xlim=c(0,1000), colour=cond, binwidth=100,xlab='Milliseconds')

where num_presentations would be 20 for condition 'a' and 10 for condition 'b' and 1000/binwidth would just be a constant to get the unit correct?

Brian Diggs · Answer 1 · 2011-08-16T21:49:57.387

2

It does not average over the conditions; it sums over them. Since condition a has 20x40 = 800 points and condition b has 10*80 = 800 points, the "area" under these "histograms" will be the same. You want each trial within a condition to have equal weight, not each point to have equal weight. That will have to be done as a pre-processing step.

trial.group <- unique(AB[,c("trial","cond")])
hists <- dlply(AB, .(trial), function(x) {hist(x$spiketime, breaks=10, plot=FALSE)})
hists.avg <- ddply(trial.group, .(cond), function(x) {
  hist.group <- ldply(hists[x$trial], function(y) {
    data.frame(mids=y$mids, counts=y$counts)
  })
  ddply(hist.group, .(mids), summarise, counts=mean(counts))
})

ggplot(data=hists.avg, aes(x=mids, y=counts, colour=cond)) + geom_line()

This is using hist to bin the data for each trial separately, then averaging the counts over the trial groups. This gives each condition equal weight, and each trial equal weight within each condition.

EDIT 1:

Taking @kohske solution, but calculating the number of trials instead of entering them explicitly:

tmp <- as.data.frame(table(unique(AB[,c("trial","cond")])["cond"]))
names(tmp) <- c("cond","ntrial")
AB <- merge(AB, tmp)

ggplot(AB, aes(spiketime, ntrial=ntrial, colour=cond)) + 
  stat_bin(aes(y=..count../ntrial*1000/100), binwidth=100, geom="line", position="identity") +
  xlim(0,1000) + 
  labs(x='Milliseconds', y="Firing rate [times/s]")

edited Aug 16 '11 at 21:49

answered Aug 15 '11 at 21:55

Brian Diggs

57,757
13
166
188

I'm sorry, there was a mistake in question! I want it to average over presentations, not conditions. Additional presentation of a condition would not vary the firing rate (which I want to plot in the histogram), it should just make the measurement less noisy. – Jonathan Aug 15 '11 at 22:15
The way you structured your data, there are always exactly 40 firings per second in condition a, and 80 firings per second in condition b. You are just varying where in the second the firings are. Are you wanting to get the rate (number of firings) within a 1/10 of a second (100 ms), and look at the average of those? If so, my solution does that. The graph for the 20 firings/sec is less noisy than the 40 (but also a lower rate; about 4 firings per 100 ms, as expected). – Brian Diggs Aug 15 '11 at 22:36
Also, I should have said averaging OVER presentations WITHIN conditions, not over conditions. – Brian Diggs Aug 15 '11 at 22:38
I was trying to avoid solving this solution through pre-processing, but it does look like it solves the problem. – Jonathan Aug 16 '11 at 20:31
I like your 'EDIT 1'. I was working on doing this myself and came up with a less elegant solution. – Jonathan Aug 16 '11 at 23:27

score 2 · Accepted Answer · answered Aug 16 '11 at 00:00

2

here is a solusion:

AB$ntrial <- ifelse(AB$cond=="a", 20, 10)
ggplot(AB, aes(spiketime, ntrial=ntrial, colour=cond)) + 
  stat_bin(aes(y=..count../ntrial*1000/100), binwidth=100, geom="line", position="identity") +
  xlim(0,1000) + 
  labs(x='Milliseconds', y="Firing rate [times/s]")

answered Aug 16 '11 at 00:00

kohske

65,572
8
165
155

Thanks a lot! This was the type of solution that I was looking for. – Jonathan Aug 16 '11 at 20:12
Unfortunately, Hadley (the author of the package) states that using extra aesthetics (ntrial) is undefined behavior. It can result in "Error in eval(expr, envir, enclos) : object 'ntrial' not found" in certain cases. – Jonathan Aug 18 '11 at 15:47

How to plot a peristimulus time histogram (PSTH) in R with ggplot2

2 Answers2