estimateSpikes {LZeroSpikeInference}

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

Description

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

Usage

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

Arguments

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)

Details

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.

Value

Returns a list with elements:

changePts the set of changepoints

spikes the set of spikes

fittedValues estimated calcium concentration

See Also

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

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

Simulation: simulateAR1, plot.simdata.

Examples

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)