# plot.transformedTrain: Plot Diagnostics for an transformedTrain Object In STAR: Spike Train Analysis with R

## Description

Six plots (selectable by `which`) are currently available: the first 5 of which correspond to Fig. 9 to 13 of Ogata (1988). The sixth one is new (as far as I know) and is still "experimental". They are all testing the first argument of `plot.transformedTrain` against the Poisson process hypothesis..

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10``` ```## S3 method for class 'transformedTrain' plot(x, which = 1:5, main, caption = c("Uniform on Trans. Obs. Time Test", "Berman's Test", "Log Survivor Function", "Lag 1 Transformed Intervals", "Variance vs Mean", "Martingale vs Trans. Time"), ask = TRUE, ...) ```

## Arguments

 `x` a `transformedTrain` object. `which` if a subset of the plots is required, specify a subset of the numbers `1:6`. `main` title to appear above the plots, if missing the corresponding element of `caption` will be used. `caption` Default caption to appear above the plots or, if `main` is given, bellow it `ask` logical; if `TRUE`, the user is asked before each plot, see `par(ask=.)`. `...` not used only there for compatibility with `plot` generic method.

## Details

If the `transformedTrain` object `x` is a the realization of a homogeneous Poisson process then, conditioned on the number of events observed, the location of the events is uniform on the (time transformed) period of observation. This is a basic property of the homogeneous Poisson process derived in Chap. 2 of Cox and Lewis (1966) and Daley and Vere-Jones (2003). This is what the first plot generated (by default) tests with a Kolmogorov-Smirnov Test. The two dotted lines on both sides of the diagonal correspond to 95 and 99% confidence intervals. This is the plot shown on Fig. 9 (p 19) of Ogata (1988).

If we write x[i] the elements of the `transformedTrain` object `x` and if the latter is the realization of a homogeneous Poisson process then the intervals:

y[i]=x[i+1]-x[i]

are iid rv from an exponential distribution with rate 1 and the:

u[i]=1 - exp(-y[i])

are iid rv from a uniform distribution on [0,1). The second plot generated (by default) tests this uniform distribution hypotheses with a Kolmogorov-Smirnov Test. This is the plot shown on Fig. 10 (p 19) of Ogata (1988) which was suggested by Berman. This is also the plot proposed by Brown et al (2002). The two dotted lines on both sides of the diagonal correspond to 95 and 99% confidence intervals.

Following the line of the previous paragraph, if the distribution of the y[i] is an exponential distribution with rate 1, then their survivor function is: exp(-y). This is what's shown on the third plot generated (by default) using a log scale for the ordinate. The point wise CI at 95 and 99% are also drawn (dotted lines). This is the plot shown on Fig. 12 (p 20) of Ogata (1988)

If the u[i] of the second paragraph are iid uniform rv on [0,1) then a plot of u[i+1] vs u[i] should fill uniformly the unit square [0,1) x [0,1). This is the fourth generated plot (by default). This is the plot shown on Fig. 11 (p 20) of Ogata (1988)

If the x[i] are realization of a homogeneous Poisson process observed between 0 and T (on the transformed time scale), then the number of events observed on non-overlapping windows of length t should be iid Poisson rv with mean t (and variance t). The observation period is chopped into non-overlapping windows of increasing length and the empirical variance of the event count is plotted versus the empirical mean, together with 95 and 99% CI (using a normal approximation). This is done by calling internally `varianceTime`. That's what's generated by the fifth plot (by default). This is the plot shown on Fig. 13 (p 20) of Ogata (1988)

The last plot is experimental and irrelevant for spike trains transformed after a `gam` or a `glm` fit. It should be useful for parametric models fitted with the maximum likelihood method.

## Author(s)

Christophe Pouzat [email protected]

## References

Cox, D. R. and Lewis, P. A. W. (1966) The Statistical Analysis of Series of Events. John Wiley and Sons.

Daley, D. J. and Vere-Jones D. (2003) An Introduction to the Theory of Point Processes. Vol. 1. Springer.

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.

`transformedTrain`, `summary.transformedTrain`, `mkGLMdf`
 ``` 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) ```