| cv.estimateSpikes {LZeroSpikeInference} | R Documentation |
Cross-validate and optimize model parameters
Description
Cross-validate and optimize model parameters
Usage
cv.estimateSpikes(dat, type = "ar1", gam = NULL, lambdas = NULL,
nLambdas = 10, hardThreshold = TRUE)
Arguments
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) |
Details
We perform cross-validation over a one-dimensional grid of \lambda values.
For each value of \lambda 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 \gamma. Given the resulting set of changepoints, we solve a constrained optimization problem for \gamma. We then refit the optimization problem with the optimized value of \gamma,
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 \gamma 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>
Value
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
See Also
Estimate spikes:
estimateSpikes,
print.estimatedSpikes,
plot.estimatedSpikes.
Cross validation:
cv.estimateSpikes,
print.cvSpike,
plot.cvSpike.
Simulation:
simulateAR1,
plot.simdata.
Examples
# 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)