peacots-package: Periodogram peaks in correlated time series
In peacots: Periodogram Peaks in Correlated Time Series

Description Details Author(s) References Examples

Calculate the periodogram of a time series, maximum-likelihood fit an Ornstein-Uhlenbeck state space (OUSS) null model to the periodogram and evaluate the statistical significance of periodogram peaks against the OUSS null hypothesis.

Package:	peacots
Type:	Package
Version:	1.2
Date:	2015-06-13
License:	GPL-3

The OUSS is a parsimonious model for a stochastically fluctuating variable (e.g. population size) with linear stabilizing forces, subject to uncorrelated measurement errors. Periodogram peaks (putative periodicities) are evaluated against this null hypothesis. In contrast to the white noise null model (the classical null model against which cyclicity is often evaluated), the OUSS process accounts for non-zero correlations between measurements and corrects for the resulting increased power at low frequencies.

Use evaluate.pm to calculate the periodogram of a time series, fit the OUSS null model and calculate the statistical significance of the periodogram maximum.

Use plotReport to generate a simple plot of the results returned by evaluate.pm.

Use significanceOfLocalPeak to evaluate the statistical significance of a secondary peak (i.e. non-global maximum) in the periodogram.

Use runExample to run an example peacots analysis based on simulation data.

Use evaluate.pm.wn to evaluate the statistical significance of the periodogram maximum against the white noise null hypothesis. This is the classical test, included for comparison.

Use ps_ouss_asymptotic to calculate the power spectrum of a particular OUSS process.

Use ps_ouss to calculate the expected periodogram from a finite time series of a particular OUSS process.

Use generate_ouss to generate random time series of a particular OUSS process.

Stilianos Louca

Louca, S., Doebeli, M. (2015) Detecting cyclicity in ecological time series, Ecology 96: 1724–1732

# Generate a cyclic time series and analyse using peacots

# Parameters
lambda         = 1;    # inverse correlation time of OU process
cyclePeriod    = 1;
cycleAmplitude = 0.6;
times          = seq(0,20,0.25);

# Example 1

# 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,
                      startRadius=2);

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


# Example 2 (using the same time series)
# In this example we don't use low-frequency trimming
# Instead, we will focus on a particular (local) periodogram peak
# and estimate its 'local' statistical significance

# calculate periodogram and fit OUSS model
report    = evaluate.pm(times=times, signal=signal, startRadius=2);

# find the periodogram mode approximately corresponding to the frequency we are interested in
cycleMode = which(report$frequencies>=0.99/cyclePeriod)[1];

# calculate local P-value for this peak
Pvalue    = significanceOfLocalPeak(power_o     = report$power_o, 
                                    lambda      = report$lambda, 
                                    power_e     = report$power_e, 
                                    time_step   = report$time_step,
                                    series_size = length(times),
                                    Nfreq       = length(report$frequencies), 
                                    peakFreq    = report$frequencies[cycleMode], 
                                    peakPower   = report$periodogram[cycleMode]);

# print result
cat(sprintf("Local P-value = %.3g for peak at frequency=%.3g\n", 
    Pvalue, report$frequencies[cycleMode]));