LZeroSpikeInference (version 1.0.1)

cv.estimateSpikes: 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)

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)

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

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.

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>

See Also

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

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

Simulation: simulateAR1, plot.simdata.

Examples

# NOT RUN {
# 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)
# }