| Title: | Exact Spike Train Inference via L0 Optimization |
|---|---|
| Description: | An implementation of algorithms described in Jewell and Witten (2017) <arXiv:1703.08644>. |
| Authors: | Sean Jewell [aut, cre] |
| Maintainer: | Sean Jewell <[email protected]> |
| License: | GPL-3 |
| Version: | 1.0.5 |
| Built: | 2025-03-01 05:36:56 UTC |
| Source: | https://github.com/jewellsean/lzerospikeinference |
Cross-validate and optimize model parameters
cv.estimateSpikes(dat, type = "ar1", gam = NULL, lambdas = NULL, nLambdas = 10, hardThreshold = TRUE)cv.estimateSpikes(dat, type = "ar1", gam = NULL, lambdas = NULL, nLambdas = 10, hardThreshold = TRUE)
dat |
fluorescence trace (a vector) |
type |
type of model, must be one of AR(1) 'ar1', or AR(1) with intercept 'intercept' |
gam |
a scalar value for the AR(1)/AR(1) + intercept decay parameter |
lambdas |
vector of tuning parameters to use in cross-validation |
nLambdas |
number of tuning parameters to estimate the model (grid of values is automatically produced) |
hardThreshold |
boolean specifying whether the calcium concentration must be non-negative (in the AR-1 problem) |
We perform cross-validation over a one-dimensional grid of values.
For each value of in this grid, we solve the corresponding optimization problem, that is, one of
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.
on a training set using a candidate value for . Given the resulting set of changepoints, we solve a constrained optimization problem for . We then refit the optimization problem with the optimized value of ,
and then evaluate the mean squared error (MSE) on a hold-out set. Note that in the final output of the algorithm,
we take the square root of the optimal value of in order to address the fact that the cross-validation
scheme makes use of training and test sets consisting of alternately-spaced timesteps.
If there is a tuning parameter lambdaT in the path [lambdaMin, lambdaMax] that produces a fit with less than 1 spike per 10,000 timesteps the path is truncated to [lambdaMin, lambdaT] and a warning is produced.
See Algorithm 3 of Jewell and Witten (2017) <arXiv:1703.08644>
A list of values corresponding to the 2-fold cross-validation:
cvError the MSE for each tuning parameter
cvSE the SE for the MSE for each tuning parameter
lambdas tuning parameters
optimalGam matrix of (optimized) parameters, rows correspond to tuning parameters, columns correspond to optimized parameter
lambdaMin tuning parameter that gives the smallest MSE
lambda1SE 1 SE tuning parameter
indexMin the index corresponding to lambdaMin
index1SE the index corresponding to lambda1SE
Estimate spikes:
estimateSpikes,
print.estimatedSpikes,
plot.estimatedSpikes.
Cross validation:
cv.estimateSpikes,
print.cvSpike,
plot.cvSpike.
Simulation:
simulateAR1,
plot.simdata.
# Not run # sim <- simulateAR1(n = 500, gam = 0.998, poisMean = 0.009, sd = 0.05, seed = 1) # plot(sim) # AR(1) model # outAR1 <- cv.estimateSpikes(sim$fl, type = "ar1") # plot(outAR1) # print(outAR1) # fit <- estimateSpikes(sim$fl, gam = outAR1$optimalGam[outAR1$index1SE, 1], # lambda = outAR1$lambda1SE, type = "ar1") # plot(fit) # print(fit) # AR(1) + intercept model # outAR1Intercept <- cv.estimateSpikes(sim$fl, type = "intercept", # lambdas = seq(0.1, 5, length.out = 10)) # plot(outAR1Intercept) # print(outAR1Intercept) # fit <- estimateSpikes(sim$fl, gam = outAR1Intercept$optimalGam[outAR1Intercept$index1SE, 1], # lambda = outAR1Intercept$lambda1SE, type = "intercept") # plot(fit) # print(fit)# Not run # sim <- simulateAR1(n = 500, gam = 0.998, poisMean = 0.009, sd = 0.05, seed = 1) # plot(sim) # AR(1) model # outAR1 <- cv.estimateSpikes(sim$fl, type = "ar1") # plot(outAR1) # print(outAR1) # fit <- estimateSpikes(sim$fl, gam = outAR1$optimalGam[outAR1$index1SE, 1], # lambda = outAR1$lambda1SE, type = "ar1") # plot(fit) # print(fit) # AR(1) + intercept model # outAR1Intercept <- cv.estimateSpikes(sim$fl, type = "intercept", # lambdas = seq(0.1, 5, length.out = 10)) # plot(outAR1Intercept) # print(outAR1Intercept) # fit <- estimateSpikes(sim$fl, gam = outAR1Intercept$optimalGam[outAR1Intercept$index1SE, 1], # lambda = outAR1Intercept$lambda1SE, type = "intercept") # plot(fit) # print(fit)
Estimate spike train, underlying calcium concentration, and changepoints based on fluorescence trace.
estimateSpikes(dat, gam, lambda, type = "ar1", calcFittedValues = TRUE, hardThreshold = FALSE, pelt = TRUE)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)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)
This package implements an algorithm for deconvolving calcium imaging data for a single neuron in order to estimate the times at which the neuron spikes.
This package implements an algorithm for deconvolving calcium imaging data for a single neuron in order to estimate the times at which the neuron spikes. 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>
Estimate spikes:
estimateSpikes,
print.estimatedSpikes,
plot.estimatedSpikes.
Cross validation:
cv.estimateSpikes,
print.cvSpike,
plot.cvSpike.
Simulation:
simulateAR1,
plot.simdata.
# To reproduce Figure 1 of Jewell and Witten (2017) <arXiv:1703.08644> sampleData <- simulateAR1(n = 500, gam = 0.998, poisMean = 0.009, sd = 0.15, seed = 8) fit <- estimateSpikes(sampleData$fl, gam = 0.998, lambda = 8, type = "ar1") plot(fit)# To reproduce Figure 1 of Jewell and Witten (2017) <arXiv:1703.08644> sampleData <- simulateAR1(n = 500, gam = 0.998, poisMean = 0.009, sd = 0.15, seed = 8) fit <- estimateSpikes(sampleData$fl, gam = 0.998, lambda = 8, type = "ar1") plot(fit)
Plot mean squared error vs. tuning parameter from the cross-validation output
## S3 method for class 'cvSpike' plot(x, ...)## S3 method for class 'cvSpike' plot(x, ...)
x |
output from cross validation procedure |
... |
arguments to be passed to methods |
Plot the solution to an L0 segmentation problem
## S3 method for class 'estimatedSpikes' plot(x, xlims = NULL, ...)## S3 method for class 'estimatedSpikes' plot(x, xlims = NULL, ...)
x |
output from running estimatedSpikes |
xlims |
optional parameter to specify the x-axis limits |
... |
arguments to be passed to methods |
Estimate spikes:
estimateSpikes,
print.estimatedSpikes,
plot.estimatedSpikes.
Cross validation:
cv.estimateSpikes,
print.cvSpike,
plot.cvSpike.
Simulation:
simulateAR1,
plot.simdata.
Plot simulated data
## S3 method for class 'simdata' plot(x, xlims = NULL, ...)## S3 method for class 'simdata' plot(x, xlims = NULL, ...)
x |
output data from simulateAR1 |
xlims |
optional parameter to specify the x-axis limits |
... |
arguments to be passed to methods |
Plot with simulated fluorescence (dark grey circles), calcium concentration (dark green line) and spikes (dark green tick marks on x-axis)
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)sim <- simulateAR1(n = 500, gam = 0.998, poisMean = 0.009, sd = 0.05, seed = 1) plot(sim)
Print CV results
## S3 method for class 'cvSpike' print(x, ...)## S3 method for class 'cvSpike' print(x, ...)
x |
output from CV |
... |
arguments to be passed to methods |
Print estimated spikes
## S3 method for class 'estimatedSpikes' print(x, ...)## S3 method for class 'estimatedSpikes' print(x, ...)
x |
estimated spikes |
... |
arguments to be passed to methods |
Print simulated data
## S3 method for class 'simdata' print(x, ...)## S3 method for class 'simdata' print(x, ...)
x |
simulated data |
... |
arguments to be passed to methods |
Simulate fluorescence trace based on simple AR(1) generative model
simulateAR1(n, gam, poisMean, sd, seed)simulateAR1(n, gam, poisMean, sd, seed)
n |
number of timesteps |
gam |
AR(1) decay rate |
poisMean |
mean for Poisson distributed spikes |
sd |
standard deviation |
seed |
random seed |
Simulate fluorescence trace based on simple AR(1) generative model
y_t = c_t + eps, eps ~ N(0, sd)
c_t = gam * c_t-1 + s_t
s_t ~ Pois(poisMean)
spikes, fluorescence, and 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)sim <- simulateAR1(n = 500, gam = 0.998, poisMean = 0.009, sd = 0.05, seed = 1) plot(sim)