如何使用 ggplot2 在 R 中绘制周刺激时间直方图 (PSTH)

Question

假设我有两个条件，“a”和“b”。神经元在条件“a”下平均每秒激发 40 个尖峰，在条件“b”下平均每秒激发 80 个尖峰。对条件“a”的响应出现 20 次，对条件“b”的响应出现 10 次，每次出现 1000 毫秒。

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))
        }
  )
)

可以使用以下方法绘制简单的直方图：

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

但是，这不会对演示文稿进行平均，因此 y 轴上的值大致相同。我想沿 y 轴绘制尖峰率（尖峰/秒），这将表明条件“b”每秒引起大约两倍的尖峰。尖峰率不会随着演示数量的增加而增加，只是变得不那么嘈杂。有没有办法在不预处理数据帧 AB 的情况下做到这一点？

换句话说，我可以按照以下方式做一些事情：

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')

其中条件“a”的 num_presentations 为 20，条件“b”的 num_presentations 为 10，并且 1000/binwidth 只是一个常量以使单位正确？

Answer 1

这是一个解决方案：

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]")

Answer 2

它不会对条件进行平均；它总结了它们。由于条件 a 有 20x40 = 800 个点，条件 b 有 10*80 = 800 个点，因此这些“直方图”下的“面积”将相同。您希望条件内的每个试验都具有相同的权重，而不是每个点都具有相同的权重。这必须作为预处理步骤来完成。

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()

这是使用 hist 分别对每个试验的数据进行分类，然后对试验组的计数进行平均。这使得每个条件具有相同的权重，并且每个条件内的每个试验具有相同的权重。

在这里，我采用 Kohske 的解决方案，但计算试验次数而不是显式输入：

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]")