thinProcess: Simulate and Analyse Data From a Model Fitted With gss

In STAR: Spike Train Analysis with R

Description

Functions thinProcess simulates a spike train and mkSimFct returns a simulating function from a gssanova fitted model. Ogata's thinning simulation method is used. Functions maxIntensity, mkSelf and mkMappedI are utility functions for the first two. Function mkPostSimAnalysis returns a function analysing a simulated spike train. Functions mkSimFct and mkPostSimAnalysis return functions which can in principle be safely used in parallel applications, that is, they have everything they need in their closure.

Usage

thinProcess(object, m2uFctList, trueData, formerSpikes,
            intensityMax, ...)
maxIntensity(object, dfWithTime, ...)
mkSelf(m2uSelf)
mkMappedI(m2uI, lag = 1)
mkSimFct(object, m2uFctList, trueData, formerSpikes,
         intensityMax, ...)
mkPostSimAnalysis(stList, trainNumber = 1, timeWindow,
                  objects, dfFct)

Arguments

object	A ssanova or a ssanova0 object.
m2uFctList	A list of functions. There should be as many functions as there are "internal" variables in object. An internal variable is a variable whose value is changed by the occurrence of a spike, like the elapsed time since the last spike of the simulated neuron, the duration of a former inter spike interval of a given lag, etc. The names of the components (functions) of the list should be the names of the variables. Each function should correspond to the map to uniform function used before fitting the model.
m2uSelf	The map to uniform function used to transform the actual elapsed time since the last spike values before fitting the model.
m2uI	The map to uniform function used to transform the actual former ISI durations before fitting the model.
lag	The considered lag (integer > 0).
trueData	A data frame containing the "true data" of the simulated epoch. This is to ensure that "external" variables such as the elapsed time since the last spike of a functionally coupled neuron are available.
formerSpikes	A vector of previous spike times. This is to make the computation of former inter spike intervals possible in every case.
intensityMax	The value of the maximal intensity. If missing function maxIntensity is called to estimate it.
dfWithTime	A data frame with one variable named "time". The latter variable is used to obtain the bin width with which the original spike train was discretized.
stList	The list of spikeTrain objects with one of the trains partly simulated. A single (partly simulated) spikeTrain object can also be used.
trainNumber	An integer, the index of the modeled and simulated spike train in stList.
timeWindow	A numeric vector of length 1 or 2. This argument specifies the time domain over which the fits contained in argument objects was performed. It is implicitly assumed that the (partial) simulation was performed outside this time domain. When a vector of length 1 is used the fitting time domain is taken as c(0,timeWindow).
objects	A list of ssanova or ssanova0 objects. Each element of the list is a "model" with which analysis will be performed. A single ssanova or ssanova0 object can also be used.
dfFct	A function whose argument is a the same as the first argument of function mkGLMdf and which returns a data frame suitable for use as argument newdata in predict.ssanova.
...	Additional arguments passed to optim which is called internally with the "BFGS" method in function maxIntensity. In functions thinProcess and mkSimFct, additional argument passed to function maxIntensity if necessary.

Details

Function thinProcess simulates a spike train with Ogata's thinning method (Ogata, 1981). The latter method required the maximal intensity of the process to be known. If such is not the case, that is, if argument intensityMax is left missing, a proposed maximal intensity is obtained with function maxIntensity. If during the simulation an actual intensity larger than the given intensityMax occurs, the simulation is interrupted and an error message is generated.

Function maxIntensity uses the central point of the variable space as its intial guess. The "BFGS" method of optim is used to find the maximal intensity.

Function mkPostSimAnalysis uses function findGlobals in order to find among the functions called by dfFct the ones which are defined in the global environment. These functions are copied in the environment (Gentleman and Ihaka, 2000) of the function returned by mkPostSimAnalysis. If the global environment defined function called by dfFct do not call themselves over functions defined in the global environment, the returned function can be safely used as argument fun of package snow's clusterApply function.

Value

thinProcess returns a spikeTrain object.

maxIntensity returns the "proposed" maximal intensity (in Hz).

mkSelf returns a function taking two arguments: self(proposedtime,st).

mkMappedI returns a function taking two arguments: function(proposedtime,st).

mkSimFct returns a function simulating a spikeTrain object. The simulation is done with function thinProcess. The returned function takes no argument. The maximal intensity required by the thinning method is stored in the closure of the returned function.

mkPostSimAnalysis returns a function taking a spikeTrain object as its single argument. This function returns a list of lists. Each list correspond to one of the models in argument objects. Each sub list has two components: lpp (the log predictive probability) and ttt (the time transformed train, a CountingProcessSamplePath object).

Note

These functions are designed to implement a rather specific type of analysis which is exposed in the "big STAR tutorial" available at: http://sites.google.com/site/spiketrainanalysiswithr/. The exemple below shows a "complete" analysis, more details and other exemples can be found in the big tutorial.

Author(s)

Christophe Pouzat [e-mail]

References

Gentleman, R. and Ihaka, R. (2000) Lexical Scope and Statistical Computing. Journal of Computational and Graphical Statistics 9: 491-508.

Ogata, Y. (1981) On Lewis' simulation method for point processes. IEEE Transactions on Information Theory IT-29: 23-31.

See Also

gssanova, as.spikeTrain, mkGLMdf, mkCPSP,

Examples

## Not run: 
## The run times given in the sequel were measured on a laptop
## with a dual core CPU: 2x Intel Core 2 Duo CPU P9500 @ 2.53GHz
## The RAM was 4 GB large. The PC ran Ubuntu 9.04 and R-2.9.2
## compiled with a dynamically link ATLAS version of BLAS.

## Start by loading the data set into the work space.
data(e060824spont)
## Get a summary of neuron 1 spike train. 
summary(e060824spont[["neuron 1"]])
## Run an automatic analysis of the train (takes ~ 4.22 s)
reportHTML(e060824spont[["neuron 1"]],filename="e060824spont_1",otherST=e060824spont[c(2)],maxiter=100)
## The renewal tests show that the discharge is not the one of
## a renewal process. The cross-correlogram shows no sign of
## coupling between the two neurons of the data set.

## Compute the partial autocorrelation function to get an idea
## of how many previous interspike intervals (ISIs) should be
## included in the model.
acf(diff(e060824spont[["neuron 1"]]),type="partial")
## The pacf plot suggests that the last ISI should be enough.
## Build the data frame.
DFA <- mkGLMdf(e060824spont[["neuron 1"]],0.004,0,59)
DFA <- within(DFA,i1 <- isi(DFA,lag=1))
DFA <- DFA[complete.cases(DFA),]

## look a the ECDF of the elapsed time since the last spike,
## that is, variable lN.1 of the data frame and of the last
## ISI (variable i1 of the data frame).
layout(matrix(1:2,nc=2))
with(DFA,plot(ecdf(lN.1),pch="."))
with(DFA,lines(range(lN.1),c(0,1),col=2,lty=2))
with(DFA,plot(ecdf(i1),pch="."))
with(DFA,lines(range(i1),c(0,1),col=2,lty=2))
## The distributions of these varaibles are clearly (and not
## surprisingly) non-uniform.

## Build emprirical functions mapping the two variables to uniform
## ones
m2u1 <- mkM2U(DFA,"lN.1",0,28.5)
m2ui <- mkM2U(DFA,"i1",0,28.5,maxiter=200)
DFA <- within(DFA,e1t <- m2u1(lN.1))
DFA <- within(DFA,i1t <- m2ui(i1))
## Cehck that the transformations did their job.
with(DFA,plot(ecdf(e1t),pch="."))
with(DFA,lines(range(e1t),c(0,1),col=2,lty=2))
with(DFA,plot(ecdf(i1t),pch="."))
with(DFA,lines(range(i1t),c(0,1),col=2,lty=2))

## The heavy computations to follow will be performed
## in parallel using the snow package.
library(snow)
## Define the number of slaves
nbSlaves <- 2
## Create the cluster.
cl <- makeCluster(rep("localhost",nbSlaves),type="SOCK")
## load STAR on each slave.
clusterEvalQ(cl,library(STAR))

## Define a function making a function performing the 
## fit with gssanova and suitable for a parallel implementation.
## The returned function does in addition time transform the
## spike train on the time window not used for the fit.
mkPFct <- function(df=DFA) {
  df
  PFct <- function(gtime,
                   fmla=event~e1t*i1t,
                   seed=20061001) {
    GF <- gssanova(fmla,
                   data=subset(df, time %in% gtime),
                   family="binomial",
                   seed=seed)
    tt <- GF %tt% subset(df, !(time %in% gtime))
    list(GF=GF,tt=tt)
  }
  PFct
}

PFct1 <- mkPFct()
## Now PFct1 returns a list with two elements: the "fit object" (GF)
## and the time transformed train (tt)

## Create a list suitable as the second argument for clusterApply
gtList <- list(early=with(DFA,time[time<=29.5]),
               late=with(DFA,time[time>29.5])
               )

## Fit and test a model with interaction between the (mapped)
## ellasped time since the last spike and the (mapped) last
## ISI. This takes ~ 95 s.
GF1.e060824spont.1 <- clusterApply(cl, gtList, PFct1)

## Look a the test battery
plot(summary(GF1.e060824spont.1[[1]][[2]]),which=c(1,2,4,6))
plot(summary(GF1.e060824spont.1[[2]][[2]]),which=c(1,2,4,6))

## Fit and test a model without interaction between the (mapped)
## ellasped time since the last spike and the (mapped) last
## ISI. This takes ~ 61 s.
GF2.e060824spont.1 <- clusterApply(cl, gtList, PFct1,fmla=event ~ e1t+i1t)

## Look a the test battery
plot(summary(GF2.e060824spont.1[[1]][[2]]),which=c(1,2,4,6))
plot(summary(GF2.e060824spont.1[[2]][[2]]),which=c(1,2,4,6))

## Compute the "predictive log probability" with Model 2
## (without interaction). This takes ~ 1.6 s
(GF2.e060824spont.1.logProb <- predictLogProb(GF2.e060824spont.1[[1]][[1]],subset(DFA,time>29.5))+predictLogProb(GF2.e060824spont.1[[2]][[1]],subset(DFA,time<=29.5)))

## Compute the "predictive log probability" with Model 1
## (with interaction). This takes ~ 3.5 s
(GF1.e060824spont.1.logProb <- predictLogProb(GF1.e060824spont.1[[1]][[1]],subset(DFA,time>29.5))+predictLogProb(GF1.e060824spont.1[[2]][[1]],subset(DFA,time<=29.5)))

## Prepare the simulations using Model 1 and 2
## Define a function initializing a mrg32k3a RNG from
## the rstream package on each slave
initMRG32k3a <- function(cl) {
  clusterEvalQ(cl,library(rstream))
  invisible(clusterCall(cl,
                        function() {
                          cmd <- parse(text=".s <- new(\"rstream.mrg32k3a\")")
                          eval(cmd,env=globalenv())
                        }
                        )
            )
  cat(paste(paste(clusterEvalQ(cl,rstream.sample(.s)),collapse=","),"\n"))
  invisible(clusterEvalQ(cl,rstream.reset(.s)))
}

## Define a function returning a list of independent and packed
## mrg32k3a rngs.
mkLecuyerList <- function(cl, ## a snow cluster
                          seed,
                          ...) {
  nbWorkers <- length(cl)
  lecuyerList <- vector(mode="list",length=nbWorkers)
  for (wIdx in 1:nbWorkers) {
    if (wIdx == 1) {
      if (!missing(seed)) lecuyerList[[1]] <- new("rstream.mrg32k3a",seed=seed)
      else lecuyerList[[1]] <- new("rstream.mrg32k3a")
    } else lecuyerList[[wIdx]] <- new("rstream.mrg32k3a")
    rstream.packed(lecuyerList[[wIdx]]) <- TRUE
  }
  lecuyerList
}

## Define a function setting the uniform rng of each slave
## to one of the independent mrg32k3a rngs created by
## mkLecuyerList.
clusterSetupRSTREAM <- function(cl, ## a snow cluster
                                lecuyerList
                                ) {


  setLecuyer <- function(packedlecuyer) {
    assign("lecuyer",packedlecuyer,env=globalenv())
    cmd <- parse(text="rstream.packed(lecuyer)<-FALSE")
    eval(cmd,env=globalenv())
  }

  clusterApply(cl,lecuyerList,setLecuyer)
  clusterEvalQ(cl,rstream.RNG(lecuyer))

}

## Load package rstream on master.
library(rstream)
## Initialize mrg32k3a rngs on each slave.
initMRG32k3a(cl)
## Create the list of independent mrg32k3a rngs
## on master.
theList <- mkLecuyerList(cl,seed=rep(20061001,6))
## Set the uniform rng of each slave to one of the
## independent mrg32k3a rngs just created.
clusterSetupRSTREAM(cl,theList)

## Define a list of map to uniform functions
fList.e060824spont.1 <- list(e1t=mkSelf(m2u1), i1t=mkMappedI(m2ui))
## Define a simulating function from Model 1 fitted on the
## half of the data set.
simF1.e060824spont.1 <- mkSimFct(object=GF1.e060824spont.1[[1]][[1]],
                                 m2uFctList=fList.e060824spont.1,
                                 trueData=subset(DFA,time>29.5),
                                 formerSpikes=with(DFA,time[event==1][time[event==1] <= 29.5])
                                 )

## Define a simulating function from Model 2 fitted on the
## half of the data set.
simF2.e060824spont.1 <- mkSimFct(object=GF2.e060824spont.1[[1]][[1]],
                                 m2uFctList=fList.e060824spont.1,
                                 trueData=subset(DFA,time>29.5),
                                 formerSpikes=with(DFA,time[event==1][time[event==1] <= 29.5])
                                 )

## Define the number of simulations to carry out.
nbRep <- 100
## Simulate spike trains in parallel using Model 1.
## This takes ~ 670 s.
sim1.e060824spont.1 <- clusterApply(cl,
                                    rep(nbRep/nbSlaves,nbSlaves),
                                    function(n,SF) lapply(1:n, function(idx) SF()),
                                    SF=simF1.e060824spont.1)

## Convert the returned list of lists into a single
## big list.
sim1.e060824spont.1 <- c(sim1.e060824spont.1[[1]],
                         sim1.e060824spont.1[[2]])

## Simulate spike trains in parallel using Model 1.
## This takes ~ 425 s.
sim2.e060824spont.1 <- clusterApply(cl,
                                    rep(nbRep/nbSlaves,nbSlaves),
                                    function(n,SF) lapply(1:n, function(idx) SF()),
                                    SF=simF2.e060824spont.1)

## Convert the returned list of lists into a single
## big list.
sim2.e060824spont.1 <- c(sim2.e060824spont.1[[1]],
                         sim2.e060824spont.1[[2]])

## Define a function generating automatically the
## proper data frame from the simulated data.
mkDF.e060824spont.1 <- function(stList) {
  DF <- mkGLMdf(stList,0.004,0,59)
  DF <- within(DF,i1 <- isi(DF,lag=1))
  DF <- DF[complete.cases(DF),]
  DF <- within(DF,e1t <- m2u1(lN.1))
  DF <- within(DF,i1t <- m2ui(i1))
  DF
}

## Define a function analysis the simulated trains with
## both Model 1 and 2.
PSAFct <- mkPostSimAnalysis(e060824spont[[1]],1,29.5,list(GF1.e060824spont.1[[1]][[1]],GF2.e060824spont.1[[1]][[1]]),mkDF.e060824spont.1)

## Analyze the simulations done with Model 1.
## This takes ~ 400 s
sim1.e060824spont.1.psa <- clusterApply(cl,sim1.e060824spont.1,PSAFct)

## Analyze the simulations done with Model 2.
## This takes ~ 400 s
sim2.e060824spont.1.psa <- clusterApply(cl,sim2.e060824spont.1,PSAFct)

## Get the log predictive probability assuming Model 1 for
## simulations done with Model 1.
sim1.e060824spont.1.lpp1 <- sapply(sim1.e060824spont.1.psa, function(l) l[[1]]$lpp)
## Get the log predictive probability assuming Model 2 for
## simulations done with Model 1.
sim1.e060824spont.1.lpp2 <- sapply(sim1.e060824spont.1.psa, function(l) l[[2]]$lpp)
## Get the log predictive probability assuming Model 1 for
## simulations done with Model 2.
sim2.e060824spont.1.lpp1 <- sapply(sim2.e060824spont.1.psa, function(l) l[[1]]$lpp)
## Get the log predictive probability assuming Model 2 for
## simulations done with Model 2.
sim2.e060824spont.1.lpp2 <- sapply(sim2.e060824spont.1.psa, function(l) l[[2]]$lpp)

## Get the observed log predictive probability with each model.
e060824spont.1.lpp1 <- predictLogProb(GF1.e060824spont.1[[1]][[1]],subset(DFA,time>29.5))
e060824spont.1.lpp2 <- predictLogProb(GF2.e060824spont.1[[1]][[1]],subset(DFA,time>29.5))

## Get the difference of observed log predictive probabilities.
e060824spont.1.lppDiff <- e060824spont.1.lpp1 - e060824spont.1.lpp2

## Look at the correlation between the log predictive probabilities
## obtained with Model 1 and 2 with data simulated with Model 1.
plot(sim1.e060824spont.1.lpp1,sim1.e060824spont.1.lpp2,main="log prob with M2 vs log prob with M1 when M1 is true",xlab="log prob with M1",ylab="log prob with M2")

## Plot the ECDF of the log predictive probabilities obtained
## with Model 1 with data simulated with Model 1.
plot(ecdf(sim1.e060824spont.1.lpp1),pch=".",main="log prob with Model 1 when Model 1 is true")
## Show the observed value of this statistic.
segments(e060824spont.1.lpp1,0,e060824spont.1.lpp1,sum(sim1.e060824spont.1.lpp1 <= e060824spont.1.lpp1)/nbRep,col=2,lwd=2)
segments(-1600,sum(sim1.e060824spont.1.lpp1 <= e060824spont.1.lpp1)/nbRep,e060824spont.1.lpp1,sum(sim1.e060824spont.1.lpp1 <= e060824spont.1.lpp1)/nbRep,col=2,lwd=2)
## Plot the ECDF of the log predictive probabilities obtained
## with Model 2 with data simulated with Model 1.
plot(ecdf(sim1.e060824spont.1.lpp2),pch=".",main="log prob with Model 2 when Model 1 is true")
## Show the observed value of this statistic.
segments(e060824spont.1.lpp2,0,e060824spont.1.lpp2,sum(sim1.e060824spont.1.lpp2 <= e060824spont.1.lpp2)/nbRep,col=2,lwd=2)
segments(-1800,sum(sim1.e060824spont.1.lpp2 <= e060824spont.1.lpp2)/nbRep,e060824spont.1.lpp2,sum(sim1.e060824spont.1.lpp2 <= e060824spont.1.lpp2)/nbRep,col=2,lwd=2)

## Plot the ECDF of the difference of the log predictive probabilities 
## obtained with data simulated with Model 1.
plot(ecdf(sim1.e060824spont.1.lpp1-sim1.e060824spont.1.lpp2),pch=".",main="log prob with M1 - log prob with M2 when M1 is true")
## Show the observed value of this statistic.
segments(e060824spont.1.lppDiff,0,e060824spont.1.lppDiff,sum(sim1.e060824spont.1.lpp1-sim1.e060824spont.1.lpp2<=e060824spont.1.lppDiff)/nbRep,col=2,lwd=2)
segments(-10,sum(sim1.e060824spont.1.lpp1-sim1.e060824spont.1.lpp2<=e060824spont.1.lppDiff)/nbRep,e060824spont.1.lppDiff,sum(sim1.e060824spont.1.lpp1-sim1.e060824spont.1.lpp2<=e060824spont.1.lppDiff)/nbRep,col=2,lwd=2)

## Look at the correlation between the log predictive probabilities
## obtained with Model 1 and 2 with data simulated with Model 2.
plot(sim2.e060824spont.1.lpp1,sim2.e060824spont.1.lpp2,main="log prob with M2 vs log prob with M1 when M2 is true",xlab="log prob with M1",ylab="log prob with M2")
## Plot the ECDF of the log predictive probabilities obtained
## with Model 1 with data simulated with Model 2.
plot(ecdf(sim2.e060824spont.1.lpp1),pch=".",main="log prob with Model 1 when Model 2 is true")
## Show the observed value of this statistic.
segments(e060824spont.1.lpp1,0,e060824spont.1.lpp1,sum(sim2.e060824spont.1.lpp1 <= e060824spont.1.lpp1)/nbRep,col=2,lwd=2)
segments(-2000,sum(sim2.e060824spont.1.lpp1 <= e060824spont.1.lpp1)/nbRep,e060824spont.1.lpp1,sum(sim2.e060824spont.1.lpp1 <= e060824spont.1.lpp1)/nbRep,col=2,lwd=2)
## Plot the ECDF of the log predictive probabilities obtained
## with Model 2 with data simulated with Model 2.
plot(ecdf(sim2.e060824spont.1.lpp2),pch=".",main="log prob with Model 2 when Model 2 is true")
## Show the observed value of this statistic.
segments(e060824spont.1.lpp2,0,e060824spont.1.lpp2,sum(sim2.e060824spont.1.lpp2 <= e060824spont.1.lpp2)/nbRep,col=2,lwd=2)
segments(-2000,sum(sim2.e060824spont.1.lpp2 <= e060824spont.1.lpp2)/nbRep,e060824spont.1.lpp2,sum(sim2.e060824spont.1.lpp2 <= e060824spont.1.lpp2)/nbRep,col=2,lwd=2)
## Plot the ECDF of the difference of the log predictive probabilities 
## obtained with data simulated with Model 1.
## Make sure that the scale of the x axis is right.
xlim=c(min(c(-e060824spont.1.lppDiff,sim2.e060824spont.1.lpp2-sim2.e060824spont.1.lpp1)),max(sim2.e060824spont.1.lpp2-sim2.e060824spont.1.lpp1))
plot(ecdf(sim2.e060824spont.1.lpp2-sim2.e060824spont.1.lpp1),pch=".",main="log prob with M2 - log prob with M1 when M2 is true",xlim=xlim)
## Show the observed value of this statistic.
points(-e060824spont.1.lppDiff,0,pch=16,col=2)

## Stop the snow cluster.
stopCluster(cl)

## End(Not run)

STAR documentation built on May 2, 2019, 11:44 a.m.