dapp: Dynamic Admixture of Poisson Process
In neuromplex: Neural Multiplexing Analysis

Description Usage Arguments Value References See Also Examples

View source: R/dynamic_neural_model-v6.R

Fits the DAPP model to binned spiking data

 
dapp(spike.counts, lengthScale = NULL, lsPrior = NULL,
     hyper = list(prec = c(1,1), sig0 = 1.87, w=c(1,1)),
     burnIn = 1e3, nsamp = 1e3, thin = 4,
     verbose = TRUE, remove.zeros = FALSE)

`spike.counts`	A list with the following items. 'Acounts': binned spike counts under condition A presented as a matrix. Rows are bins, columns are replicates (trials). 'Bcount': binned spike counts under condition B. 'ABcounts': binned spike counts under condition AB. 'bin.mids': an array giving the mid-points of the time bins. 'bin.width': a scalar giving the bin width.
`lengthScale`	an array giving the length scale parameter values to be used for Gaussian process prior. Defaults to `sort(0.16 * resp.horiz / c(4, 3, 2, 1, 0.5, 0.1))` where `resp.horiz` is the time horizon of the response period.
`lsPrior`	an array of the same length as `lengthScale` giving the prior probabilities of the length scale values.
`hyper`	a list of hyper parameters with the following iterms. 'prec': a 2-vector giving the shape and rate parameters of the gamma distribution on the Dirichlet precision parameter. 'sig0': a scalaer giving the scale of the (centered) logistic distribution used in transforming the Gaussian random curves into curves restricted between 0 and 1.
`burnIn`	number of MCMC iterations to discard as burn-in.
`nsamp`	number of MCMC draws to be saved for posterior inference.
`thin`	the thinning rate at which MCMC draws are to be saved. The total number of iterations equals `burnIn + nsamp * thin`
`verbose`	logical indicating if some fit details should be printed during the course of the MCMC
`remove.zeros`	logical indicating if trials with zero spike count shuold be removed from the analysis

Returns a list of class "dapp" containting the following items.

`lsProb`	posterior preditctive draws of length scale
`lambda.A`	posterior draws of lambda.A at bin mid-points
`lambda.B`	posterior draws of lambda.B at bin mid-points
`alpha`	posterior draws of the alpha curves at bin mid-points
`A`	posterior draws of the latent variable A which gives the AB spike counts (by bin) that are to be attributed to signal A (the remaining are attributed to signal B)
`prec`	posterior draws of precision
`alpha.pred`	posterior predictive draws of alpha (of a future trial)
`psl.pred`	posterior predictive draw of the feature parameters (phi, psi, ell) (of a future trial)
`details`	mcmc details given as an array of `c(niter, nsamp, burnIn, thin, MH acceptance rate)`
`hyper`	hyper parameters used in model fitting
`lengthScale`	length scale set used in model fitting
`lsPrior`	length scale prior
`bin.mids`	bin mid-points
`bin.width`	bin width
`mcmc`	mcmc controls (burn-in length, thinning rate and number of saved draws)

Glynn, C., Tokdar, S.T., Zaman, A., Caruso, V.C., Mohl, J.T., Willett, S.M., and Groh, J.M. (2020+). Analyzing second order stochasticity of neural spiking under stimuli-bundle exposure. The Annals of Applied Statistics. Accepted.

plot.dapp, summary.dapp and predict.dapp.

  ## Note:
  #### The example below uses a simpler synthetic data, a wider bin-width
  #### and a shorter MCMC run to keep the run length less than 5s
  #### Use ?plot.dapp or ?plot.summary for a more realistic example
  
  ## Generate 30 A and 30 B trials with rate functions
  ##    lambda.A(t) = 160*exp(-2*t/1000) + 40*exp(-0.2*t/1000)
  ##    lambda.B(t) = 40*exp(-2*t/1000)
  ## where time t is measured in ms. Then, generate 25 AB trials,
  ## roughly 2/3 with flat weight curves with a constant intensity
  ## either close to A, or close to B (equally likely). The 
  ## remaining 1/3 curves are sinusoidal that snake between 0.01 and 0.99 
  ## with a period randomly drawn between 400 and 1000
  
  ntrials <- c(nA=30, nB=30, nAB=25)
  flat.range <- list(A=c(0.85, 0.95),
                     B=c(0.05, 0.15))
  flat.mix <- c(A=1/2, B=1/2)
  wavy.span <- c(0.01, 0.99)
  wavy.period <- c(400, 1000)
  
  T.horiz <- 1000
  rateB <- 40 * exp(-2*(1:T.horiz)/T.horiz)
  rateA <- 4*rateB + 40 * exp(-0.2*(1:T.horiz)/T.horiz)
  
  synth.data <- synthesis.dapp(ntrials = ntrials, pr.flat = 2/3,
                               intervals = flat.range, wts = flat.mix,
                               span = wavy.span, period.range = wavy.period,
                               lambda.A=rateA, lambda.B=rateB)
  
  ## Generate binned spike counts witb 100 ms bins
  spike.counts <- mplex.preprocess(synth.data$spiketimes, bw=100, visualize=FALSE)
  
  ## Fit the DAPP model to data
  #### A short MCMC run is done below to keep the run length short.
  #### Use default or larger values for burn, nsamp and thin
  #### for more reliable estimation
  fit.post <- dapp(spike.counts, burn=10, nsamp=90, thin=1, verbose=FALSE)
  
  ## Visualize model fit
  plot(fit.post)
  
  ## Post process results to assign second order stochasticity labels
  summary(fit.post)