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

Question

假设我有两个条件，“a”和“b”。一个神经元在“a”条件下平均每秒发射 40 个脉冲（Hz），在“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')

其中 num_presentations 对于条件“a”为 20，对于条件“b”为 10，而 1000/binwidth 将只是一个常数以使单位正确？

Answer 1

它不会在条件下进行平均；它总结了他们。由于条件 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 分别对每个试验的数据进行分箱，然后对试验组的计数进行平均。这使每个条件具有相同的权重，并且每个试验在每个条件内具有相同的权重。

编辑 1：

采用@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]")

Answer 2

这是一个解决方案：

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