peacots (version 1.0)

evaluate.pm: Evaluate the periodogram maximum of a time series

Description

Calculate the Lomb-Scargle periodogram of a time series and estimate the statistical significance of the periodogram maximum based on the null hypothesis of an Ornstein-Uhlenbeck state space (OUSS) process.

Usage

evaluate.pm(times, signal, minPeakFreq = 0, minFitFreq = 0, 
            iterations = 100, accuracy = 0.002, startRadius = 2, 
            verbose = FALSE)

Arguments

times

Numerical vector. Time points of the time series.

signal

Numerical vector of same size as times. Values of the time series.

minPeakFreq

Non-negative number. Minimum considered frequency when determining periodogram peak. Use this if you want to ignore low-frequency components in the spectrum when searching for the periodogram maximum.

minFitFreq

Non-negative number. Minimum considered frequency when fitting the OUSS model to the periodogram. Use this to ignore low-frequency components in the spectrum when estimating the OUSS parameters, e.g. if you suspect your time series to be trended. mi

iterations

Number of iterations for the maximum-likelihood fitter. Increasing this will result in greater estimation accuracy but also reduced performance.

startRadius

Single integer. The number of initial guesses for each OUSS parameter during optimization of the likelihood function, per direction (up and down). Increasing this will improve estimation accuracy. However, execution time scales with startRadius^2

accuracy

An upper bound for the standard deviation of the P-value estimator using repeated Bernoulli trials. The number of trials scales as 1/accuracy^2. If you are using 0.05 as a nominal significance threshold, then accuracy=0.002

verbose

Single logical. If TRUE, the function prints out occasional progress reports.

Value

• A list with the following entries:
• error
  Will be TRUE if an error occured, FALSE otherwise.
• errorMessage
  A short error message if error==TRUE.
• If error==FALSE, the returned list also includes:
• frequencies
  Available periodogram frequencies as a numerical vector.
• periodogram
  Periodogram powers corresponding to the returned frequencies, as a numerical vector.
• fittedPS
  Maximum-likelihood fitted OUSS power spectrum corresponding to frequencies, as a numerical vector.
• power_o
  Estimated power at zero-frequency generated by the underlying OU process.
• lambda
  Estimated resilience of the OU process (in inverse time units).
• power_e
  Estimated power at large frequencies due to the random measurement errors.
• time_step
  The average time step of the time series, as used to fit the OUSS power spectrum.
• peakMode
  An integer indicating the position of the periodogram maximum in the vector frequencies.
• minPeakMode
  The minimum periodogram mode considered for determining the periodogram maximum. This will be 1 if minPeakFreq==0.
• minFitMode
  The minimum periodogram mode considered for estimating the white noise power. This will be 1 if minFitFreq==0.
• MLL
  Log-likelihood value at calculated maximum.
• P
  Statistical significance of the periodogram peak against the null hypothesis of the estimated OUSS process. This is the probability that the estimated OUSS process would generate a periodogram with global maximum at least as ``extreme'' as the observed peak (among the considered frequencies). See details above.
• Plocal
  Statistical significance of the relative power of the periodogram peak. Mainly used for comparison to the statistical significance of secondary peaks. See significanceOfLocalPeak.

Details

The OUSS model describes the measurement of an Ornstein-Uhlenbeck (OU) stochastic process at discrete times with additional uncorrelated Gaussian measurement errors of fixed variance. The OU process itself is a continuous-time random walk with linear stabilizing forces, described by the stochastic differential equation $$dX = \lambda(\mu-X) dt + s dW,$$ where $W$ is the standard Wiener process. The OUSS process is a parsimonious model for describing stochastically fluctuating variables subject to linear stabilizing forces and uncorrelated measurement errors.

Due to temporal correlations, the OUSS power spectrum increases gradually towards lower frequencies, as opposed to the white noise spectrum, which is flat. Using white noise as a null hypothesis for the evaluation of cyclicity in time series, particularly for systems with long correlation times, may result in increased false cycle detection rates because of the increased low-frequency power. The OUSS model is an attempt to account for these correlations.

The OUSS model parameters are estimated using maximum-likelihood fitting to the periodogram. The likelihood function, and therefore the OUSS parameter estimates, are only approximations that become exact for long regular time series. The statistical significance of the periodogram peak is defined as the probability that the estimated OUSS process would generate a periodogram whose maximum (a) is at least as strong as the observed peak, and (b) has a power-to-expectation ratio at least as strong as the observed power-to-expectation ratio of the peak. The P-value is estimated using repeated Bernoulli trials.

If you want to evaluate secondary peaks (i.e. non-global periodogram maxima), you will need to either (a) adjust the parameters minPeakFreq and minFitFreq to omit low-frequency modes or (b) use significanceOfLocalPeak after using evaluate.pm.

References

Louca, S., Doebeli, M. - Detecting cyclicity in ecological time series (in review, as of June 2014)

See Also

evaluate.pm.wn, significanceOfLocalPeak, ps_ouss

Examples

# In this example we generate a cyclic time series
# and analyse its periodogram using evaluate.pm

# Parameters
lambda         = 1;    # inverse correlation time of OU process
cyclePeriod    = 1;    # Cycle period
cycleAmplitude = 0.75;
times          = seq(0,20,0.2);

# generate cyclic time series by adding a periodic signal to an OUSS process
signal = cycleAmplitude * cos(2*pi*times/cyclePeriod) +
         generate_ouss(times, mu=0, sigma=1, lambda=lambda, epsilon=0.5);

# Find periodogram peak and estimate statistical significance
# Ignore frequencies lower than a pre-defined threshold
# to avoid masking by low-frequency maximum
report = evaluate.pm(times=times, signal=signal, 
                     minPeakFreq=lambda/3, 
                     minFitFreq=lambda/3);

# plot overview of periodogram peak analysis
plotReport(sprintf("Cyclic at frequency %.2g",1/cyclePeriod), 
           times=times, signal=signal, report=report);