estimateSpikes: Estimate spike train, underlying calcium concentration, and...
In jewellsean/LZeroSpikeInference: Exact Spike Train Inference via L0 Optimization

Description Usage Arguments Details Value See Also Examples

Estimate spike train, underlying calcium concentration, and changepoints based on fluorescence trace.

1 2	estimateSpikes(dat, gam, lambda, type = "ar1", calcFittedValues = TRUE, hardThreshold = FALSE, pelt = TRUE)

`dat`	fluorescence data
`gam`	a scalar value for the AR(1)/AR(1) + intercept decay parameter.
`lambda`	tuning parameter lambda
`type`	type of model, must be one of AR(1) 'ar1', AR(1) + intercept 'intercept'.
`calcFittedValues`	TRUE to calculate fitted values.
`hardThreshold`	boolean specifying whether the calcium concentration must be non-negative (in the AR-1 problem)
`pelt`	boolean specifying whether PELT (default) or optimal partitioning algorithm is used to compute the segmentation. Both yield the same solution, however, PELT can be orders of magnitude faster.

This algorithm solves the optimization problems

AR(1)-model: minimize_c1,...,cT 0.5 sum_t=1^T ( y_t - c_t )^2 + lambda sum_t=2^T 1_c_t neq gamma c_t-1 for the global optimum, where $y_t$ is the observed fluorescence at the tth timepoint.

If hardThreshold = T then the additional constraint c_t >= 0 is added to the optimzation problem above.

AR(1) with intercept: minimize_c1,...,cT,b1,...,bT 0.5 sum_t=1^T (y_t - c_t - b_t)^2 + lambda sum_t=2^T 1_c_t neq gamma c_t-1, b_t neq b_t-1 where the indicator variable 1_(A,B) equals 1 if the event A cup B holds, and equals zero otherwise.

See Jewell and Witten (2017) <arXiv:1703.08644>, section 2 and 5.

Note that "changePts" and "spikes" differ by one index due to a quirk of the conventions used in the changepoint literature and the definition of a spike.

Returns a list with elements:

changePts the set of changepoints

spikes the set of spikes

fittedValues estimated calcium concentration

Estimate spikes: estimateSpikes, print.estimatedSpikes, plot.estimatedSpikes.

Cross validation: cv.estimateSpikes, print.cvSpike, plot.cvSpike.

Simulation: simulateAR1, plot.simdata.

sim <- simulateAR1(n = 500, gam = 0.998, poisMean = 0.009, sd = 0.05, seed = 1)
plot(sim)

# AR(1) model

fit <- estimateSpikes(sim$fl, gam = 0.998, lambda = 1, type = "ar1")
plot(fit)
print(fit)

# AR(1) + intercept model
fit <- estimateSpikes(sim$fl, gam = 0.998, lambda = 1, type = "intercept")
plot(fit)
print(fit)