Transform spike times from a
gam fitted model as
defined by Ogata (1988) and Brown et al (2002). If the model structure
is "correct" and if the model parameters are properly estimated the
result of the time transformation should be the realization of a
Poisson process with rate 1.
An object returned by
A binary (0,1) vector of integers with the same length
A character string defining a condition to be fulfilled
by the event in order to be selected, like:
fitted.values component of
obj contains the
(estimated) probability to observe a spike in each time bin where the covariates
required by the fitted model were defined. It is then straightforward
to show using the concept of product integral (Kalbfleisch and
Prentice, 2002; Andersen et al, 1993),provided that the time bin width
is small enough to have a very small probability in each bin, that the
cumulated sum of these probabilities is the expected number of events
observed up to a given time. This expected number of events
which is returned by
transformedTrain. It is also the result of
the "time transformation" proposed by Ogata (1988) and brought to the
spike train analysis field under the name "time rescaling (theorem)"
by Brown et al (2002).
transformedTrain can also be used to transform the times of the
spikes of neurons whose spike trains were simultaneously recorded and
discretized in exactly the same way as the neuron used to
obj. This is useful to explore the possibility of
functional interactions between a putative pre-synaptic neuron (whose
spike train would correspond to argument
target) and a
post-synaptic one used to generate
transformedTrain returns an object of class
transformedTrain inheriting from class
object is fundamentally a numeric vector with strictly increasing
elements containing the transformed times (or the expected number of events).
As mentioned only the spikes for which the covariates of the model are
available have their times transformed. That practically means that
the length of the
transformedTrain object returned by function
transformedTrain can be shorter than the length of the
spikeTrain object (or more precisely than the number
of spikes defined in
target). If one works with a model
involving the elapsed times since the last three spikes then the
fourth spike of the train will be the first to be transformed. You
should therefore expect some left truncation of the data at the
beginning of each acquisition epoch.
Christophe Pouzat email@example.com
Ogata, Yosihiko (1988) Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes. Journal of the American Statistical Association 83: 9-27.
Brown, E. N., Barbieri, R., Ventura, V., Kass, R. E. and Frank, L. M. (2002) The time-rescaling theorem and its application to neural spike train data analysis. Neural Computation 14: 325-346.
Kalbfleisch, John D. and Prentice, Ross L. (2002) The Statistical Analysis of Failure Time Data. Wiley Interscience.
Andersen, Per Kragh, Borgan, Ornulf, Gill, Richard D. and Keiding, Niels (1993) Statistical Models Based on Counting Processes. Springer-Verlag.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
## Not run: ## Let us consider neuron 1 of the CAL2S data set data(CAL2S) CAL2S <- lapply(CAL2S,as.spikeTrain) CAL2S[["neuron 1"]] renewalTestPlot(CAL2S[["neuron 1"]]) summary(CAL2S[["neuron 1"]]) ## Make a data frame with a 4 ms time resolution cal2Sdf <- mkGLMdf(CAL2S,0.004,0,60) ## keep the part relative to neuron 1, 2 and 3 separately n1.cal2sDF <- cal2Sdf[cal2Sdf$neuron=="1",] n2.cal2sDF <- cal2Sdf[cal2Sdf$neuron=="2",] n3.cal2sDF <- cal2Sdf[cal2Sdf$neuron=="3",] ## remove unnecessary data rm(cal2Sdf) ## Extract the elapsed time since the second to last and ## third to last for neuron 1. Normalise the result. n1.cal2sDF[c("rlN.1","rsN.1","rtN.1")] <- brt4df(n1.cal2sDF,"lN.1",2,c("rlN.1","rsN.1","rtN.1")) ## load mgcv library library(mgcv) ## fit a model with a tensorial product involving the last ## three spikes and using a cubic spline basis for the last two ## To gain time use a fixed df regression spline n1S.fitA <- gam(event ~ te(rlN.1,rsN.1,bs="cr",fx=TRUE) + rtN.1,data=n1.cal2sDF,family=binomial(link="logit")) ## transform time N1.Lambda <- transformedTrain(n1S.fitA) ## check out the resulting spike train using the fact ## that transformedTrain objects inherit from spikeTrain ## objects N1.Lambda ## Use more formal checks summary(N1.Lambda) plot(N1.Lambda,which=c(1,2,4,5),ask=FALSE) ## Transform spike trains of neuron 2 and 3 N2.Lambda <- transformedTrain(n1S.fitA,n2.cal2sDF$event) N3.Lambda <- transformedTrain(n1S.fitA,n3.cal2sDF$event) ## Check interactions summary(N2.Lambda %frt% N1.Lambda) summary(N3.Lambda %frt% N1.Lambda) plot(N2.Lambda %frt% N1.Lambda,ask=FALSE) plot(N3.Lambda %frt% N1.Lambda,ask=FALSE) ## End(Not run)