Как построить перистимульную гистограмму времени (PSTH) в R с помощью ggplot2 - PullRequest
Новая
       12

Как построить перистимульную гистограмму времени (PSTH) в R с помощью ggplot2

4 голосов
/ 16 августа 2011

Скажите, что у меня есть два условия: «а» и «б».Нейрон запускает в среднем 40 пиков в секунду (Гц) в состоянии «а» и 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, которая показала бы, что условие '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 будет 20 для условия'a' и 10 для условия 'b' и 1000 / binwidth будут просто константами для правильной единицы измерения?

Ответы [ 2 ]

2 голосов
/ 16 августа 2011

вот решение: 

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

2 голосов
/ 16 августа 2011

не усредняется по условиям; это их сумма. Поскольку условие 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]")
...