Package 'peacots' reference manual

Title:	Periodogram Peaks in Correlated Time Series
Description:	Calculates the periodogram of a time series, maximum-likelihood fits an Ornstein-Uhlenbeck state space (OUSS) null model and evaluates the statistical significance of periodogram peaks against the OUSS null hypothesis. The OUSS is a parsimonious model for stochastically fluctuating variables with linear stabilizing forces, subject to uncorrelated measurement errors. Contrary to the classical white noise null model for detecting cyclicity, the OUSS model can account for temporal correlations typically occurring in ecological and geological time series. Citation: Louca, Stilianos and Doebeli, Michael (2015) <doi:10.1890/14-0126.1>.
Authors:	Stilianos Louca [aut, cre, cph]
Maintainer:	Stilianos Louca <[email protected]>
License:	GPL-3
Version:	1.3.2
Built:	2024-11-13 04:58:18 UTC
Source:	https://github.com/cran/peacots

Periodogram peaks in correlated time series

Description

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.

Details

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.

Author(s)

Stilianos Louca

References

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

Examples

# 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]));
# 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]));

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.005, startRadius = 5, 
            verbose = FALSE)
evaluate.pm(times, signal, minPeakFreq = 0, minFitFreq = 0, 
            iterations = 100, accuracy = 0.005, startRadius = 5, 
            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. `minFitFreq` might differ from `minPeakFreq`.
`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.005` is probably good enough.
`verbose`	Single logical. If `TRUE`, the function prints out occasional progress reports.

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 (power $S$ at frequency $F$ ) under the null-hypothesis of an OUSS process is defined as the probability that the same OUSS process would generate a periodogram whose maximum (power $s$ at frequency $f$ ) satisfies

$s^2/e\ge S^2/E,$

where $e$ and $E$ are the expected periodogram powers at the frequencies $f$ and $F$ , respectively. The P-value is estimated via repeated Bernoulli trials in which random OUSS periodograms are emulated by exponentially distributed numbers.

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.

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 spectrum at zero-frequency generated by the underlying OU process.
`lambda`	Estimated resilience of the OU process (in inverse time units).
`power_e`	Estimated power spectrum at large frequencies due to the random measurement errors.
`sigma`	Estimated stationary standard deviation of the underlying OU process.
`epsilon`	Estimated standard deviation of 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`.

Author(s)

Stilianos Louca

References

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

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.6;
times          = seq(0,20,0.25);

# 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);     
      
# 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.6;
times          = seq(0,20,0.25);

# 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);

Statistical significance of periodogram peaks (classical)

Description

Calculates the Lomb-Scargle periodogram for the given time series and estimates the statistical significance of the global periodogram maximum based on the null hypothesis of uncorrelated (white) noise. Available for historical reasons and for comparison purposes.

Usage

evaluate.pm.wn(times, signal, minPeakFreq=0, minFitFreq=0)
evaluate.pm.wn(times, signal, minPeakFreq=0, minFitFreq=0)

Arguments

`times`	Numerical vector. Time points of the time series.
`signal`	Numerical vector of same size as `times`. Values of the time series.
`minPeakFreq`	Single non-negative number. Minimum considered frequency when determining periodogram peak. Use this to ignore low-frequency components from the spectrum.
`minFitFreq`	Single non-negative number. Minimum considered frequency when fitting the white noise null model to the periodogram. Use this to ignore low-frequency components from the spectrum.

Value

A list with the 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.
`peakMode`	An integer indicating the position of the global periodogram maximum (starting at `minPeakFreq`) in the vector `frequencies`.
`powerEstimate`	The estimated white noise power. Estimated from the average periodogram power, which corresponds to using the total variance of the time series (if `minFitFreq==0`).
`minPeakMode`	The minimum periodogram mode considered for determining the periodogram peak. 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`.
`RSS`	The sum of squared residuals of the periodogram from the estimated white noise power.
`P`	Statistical significance of periodogram peak. This is the probability that a white noise periodogram (of the estimated power) would generate a peak at least as strong as the observed peak (among the considered frequencies). The calculated P-value is only an approximation that becomes exact for long regular time series.

Author(s)

Stilianos Louca

References

Scargle, J. D. (1982) - Studies in astronomical time series analysis. II Statistical aspects of spectral analysis of unevenly spaced data, The Astrophysical Journal 263, 835–853

Horne, J. H., Baliunas, S. L. (1986) - A prescription for period analysis of unevenly sampled time series, The Astrophysical Journal 302, 757–763

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

Examples

# generate time series
times  = seq(0,20,0.25);
signal = rnorm(n=length(times));
report = evaluate.pm.wn(times=times, signal=signal);

# plot time series
old.par <- par(mfrow=c(1, 2));
plot(ts(times), ts(signal), 
     xy.label=FALSE, type="l", 
     ylab="signal", xlab="time", main="OUSS time series");
     
# plot periodogram
title = sprintf("Periodogram OUSS analysis\n(peak freq=%.3g, P=%.2g)",
                report$frequencies[report$peakMode],report$P);
plot(ts(report$frequencies),
     ts(report$periodogram), 
     xy.label=FALSE, type="l", 
     ylab="power", xlab="frequency", 
     main=title, col="black");
	
# plot fitted flat WN power
lines(c(report$frequencies[1],tail(report$frequencies,1)), 
      c(report$powerEstimate, report$powerEstimate ), col="blue");
par(old.par)

# generate time series
times  = seq(0,20,0.25);
signal = rnorm(n=length(times));
report = evaluate.pm.wn(times=times, signal=signal);

# plot time series
old.par <- par(mfrow=c(1, 2));
plot(ts(times), ts(signal), 
     xy.label=FALSE, type="l", 
     ylab="signal", xlab="time", main="OUSS time series");
     
# plot periodogram
title = sprintf("Periodogram OUSS analysis\n(peak freq=%.3g, P=%.2g)",
                report$frequencies[report$peakMode],report$P);
plot(ts(report$frequencies),
     ts(report$periodogram), 
     xy.label=FALSE, type="l", 
     ylab="power", xlab="frequency", 
     main=title, col="black");
	
# plot fitted flat WN power
lines(c(report$frequencies[1],tail(report$frequencies,1)), 
      c(report$powerEstimate, report$powerEstimate ), col="blue");
par(old.par)

Generate random time series of the OUSS process

Description

Generate a random time series of the 1-dimensional stationary Ornstein-Uhlenbeck state space (OUSS) process.

Usage

generate_ouss(times, mu, power_o, sigma, 
              lambda, power_e, epsilon)
generate_ouss(times, mu, power_o, sigma, 
              lambda, power_e, epsilon)

Arguments

`times`	Numeric vector of times for which to evaluate OUSS model. Times need to be strictly increasing.
`mu`	Single number. Deterministic equilibrium of OU process, i.e., the expected value of the time series at any particular time.
`sigma`	Single number. Standard deviation of OU fluctuations around equilibrium.
`power_o`	Single non-negative number. Power spectrum at zero-frequency generated by the OU process. Either `power_o` or `sigma` (but not both) must be provided.
`lambda`	Single non-negative number. Resilience (also known as relaxation rate) of the OU process. This is the inverse of the OU correlation time.
`epsilon`	Single number. Standard deviation of Gaussian measurement error. Setting this to zero will yield a time series from the classical OU process.
`power_e`	Single non-negative number. Asymptotic power spectrum at large frequencies due to the Gaussian measurement errors. Setting this to zero will yield a classical OU process. Either `power_e` or `epsilon` (but not both) must be provided.

Details

The OUSS model describes the measurement of an Ornstein-Uhlenbeck (OU) stochastic process at discrete times with additional uncorrelated Gaussian measurement errors. The OU process itself is a continuous-time random walk (Brownian motion) with linear stabilizing forces, described by the stochastic differential equation

$dX = \lambda(\mu-X) dt + s dW,$

where $W$ is the standard Wiener process and $s^2=2\lambda\sigma^2.$ The OUSS model is obtained by adding uncorrelated Gaussian numbers with zero mean and variance $\epsilon^2$ to the time series.

Value

A numeric vector of same length as times, containing sampled values of the OUSS process. These values will all have the same expectation (mu) and variance (sigma^2+epsilon^2) but will be correlated.

Author(s)

Stilianos Louca

References

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

Dennis, B., Ponciano, J.M. - Density dependent state-space model for population abundance data with unequal time intervals, Ecology (in press as of June 2014)

Examples

# define times
times = seq(0,100,0.5);

# generate OUSS time series
signal 	= generate_ouss(times=times, mu=0, sigma=1, lambda=1, epsilon=0.5);

# plot time series
plot(ts(times), ts(signal), 
     xy.label=FALSE, type="l", 
     ylab="signal", xlab="time", main="OUSS time series");
# define times
times = seq(0,100,0.5);

# generate OUSS time series
signal 	= generate_ouss(times=times, mu=0, sigma=1, lambda=1, epsilon=0.5);

# plot time series
plot(ts(times), ts(signal), 
     xy.label=FALSE, type="l", 
     ylab="signal", xlab="time", main="OUSS time series");

Plot results of a `evaluate.pm` analysis

Description

Create a simple plot of a time series and the results of a evaluate.pm analysis (including the periodogram and the fitted OUSS power spectrum).

Usage

plotReport(name="", times=NULL, signal=NULL, 
           report=NULL, plotFile=NULL, dataFile=NULL, sep=" ")
plotReport(name="", times=NULL, signal=NULL, 
           report=NULL, plotFile=NULL, dataFile=NULL, sep=" ")

Arguments

`name`	Character. A short name for the time series to be used for the plots (e.g. 'long-term study' or 'hare population').
`times`	Numeric vector. The time points of the time series used for the analysis. Set to `NULL` to not plot the time series.
`signal`	Numeric vector. The time series values (signal) used for the analysis. Set to `NULL` to not plot the time series.
`report`	The value returned by `evaluate.pm`. This is a list with the periodogram analysis results. Set to `NULL` to not plot any analysis results.
`plotFile`	An optional path to a PDF file to be generated with the plot.
`dataFile`	An optional path to a data file for storing the time series and the results of the analysis.
`sep`	Separator to be used for the data file. Only relevant if `dataFile` is not `NULL`.

Value

No return value.

Author(s)

Stilianos Louca

References

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

Examples

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

# find periodogram peak and estimate statistical significance
report = evaluate.pm( times=times, 
                      signal=signal, 
                      minPeakFreq=0.3, 
                      minFitFreq=0.3,
                      startRadius=2);

# plot overview of periodogram peak analysis
plotReport(sprintf("Example"), times=times, signal=signal, report=report);
# generate cyclic time series by adding a periodic signal to an OUSS process
times  = seq(0,20,0.25);
signal = 0.6 * cos(2*pi*times) + generate_ouss(times, mu=0, sigma=1, lambda=1, epsilon=0.5);

# find periodogram peak and estimate statistical significance
report = evaluate.pm( times=times, 
                      signal=signal, 
                      minPeakFreq=0.3, 
                      minFitFreq=0.3,
                      startRadius=2);

# plot overview of periodogram peak analysis
plotReport(sprintf("Example"), times=times, signal=signal, report=report);

Expected periodogram of the OUSS process

Description

Returns the expected periodogram power of the Ornstein-Uhlenbeck state space (OUSS) process at a particular frequency, when sampled at regular time intervals for a finite time.

Usage

ps_ouss(freq, power_o, sigma, 
        rho, lambda, 
        power_e, epsilon, 
        time_step, series_size)
ps_ouss(freq, power_o, sigma, 
        rho, lambda, 
        power_e, epsilon, 
        time_step, series_size)

Arguments

`freq`	Single number or numeric vector. The frequency for which to the power spectrum is to be calculated.
`power_o`	Single non-negative number. Power at zero-frequency generated by the underlying OU process, when sampled at the given `time_step`. Either `power_o` or `sigma` (but not both) must be provided.
`sigma`	Single number. Standard deviation of OU fluctuations around equilibrium. Either `power_o` or `sigma` (but not both) must be provided.
`rho`	Single number between 0 (exclusive) and 1 (inclusive). Correlation of the OU process between two subsequent time points. Either `rho` or `lambda` (but not both) must be provided.
`lambda`	Single non-negative number. Resilience (or relaxation rate) of the OU process. This is also the inverse correlation time of the OU process. Either `rho` or `lambda` (but not both) must be provided.
`power_e`	Single non-negative number. Asymptotic power at large frequencies due to the random measurement errors. Setting this to zero corresponds to the classical OU process. Either `power_e` or `epsilon` (but not both) must be provided.
`epsilon`	Single number. Standard deviation of Gaussian measurement error. Setting this to zero corresponds to the classical OU process. Either `power_e` or `epsilon` (but not both) must be provided.
`time_step`	Positive number. The time step of the time series that was (or will be) used for periodogram generation.
`series_size`	Positive integer. The number of sampled time points.

Details

The OUSS parameters power_o, lambda and power_e will typically be maximum-likelihood fitted values returned by evaluate.pm. The value of time_step is also returned by evaluate.pm and is inferred from the analysed time series. More generally, power_o and power_e are proportional to the OUSS parameters sigma^2 and epsilon^2 (see generate_ouss), respectively, but the exact scaling depends on the normalization used for the periodogram.

In the limit where series_size becomes very large, ps_ouss becomes the same as ps_ouss_asymptotic.

Value

Returns a numeric vector of the same size as freq, containing the corresponding expected periodogram powers of the OUSS process.

Author(s)

Stilianos Louca

References

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

Examples

# generate OUSS time series
times  = seq(0,20,0.25);
signal = generate_ouss(times, mu=0, sigma=1, lambda=1, epsilon=0.5);

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

# plot periodogram
plot(report$frequencies, report$periodogram, 
     type="l", ylab="power", xlab="frequency", main="periodogram & fitted OUSS power spectrum");

# plot expected OUSS periodogram
lines(report$frequencies, 
      ps_ouss(freq=report$frequencies, 
              power_o=report$power_o, 
              lambda=report$lambda, 
              power_e=report$power_e, 
              time_step=report$time_step,
              series_size=length(times)), 
      col="red");

# generate OUSS time series
times  = seq(0,20,0.25);
signal = generate_ouss(times, mu=0, sigma=1, lambda=1, epsilon=0.5);

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

# plot periodogram
plot(report$frequencies, report$periodogram, 
     type="l", ylab="power", xlab="frequency", main="periodogram & fitted OUSS power spectrum");

# plot expected OUSS periodogram
lines(report$frequencies, 
      ps_ouss(freq=report$frequencies, 
              power_o=report$power_o, 
              lambda=report$lambda, 
              power_e=report$power_e, 
              time_step=report$time_step,
              series_size=length(times)), 
      col="red");

Power spectrum of the OUSS process

Description

Returns the power spectrum of the Ornstein-Uhlenbeck state space (OUSS) process at a particular frequency. This is the asymptotic expected periodogram power for long regular time series.

Usage

ps_ouss_asymptotic(freq, power_o, sigma, 
                   rho, lambda, 
                   power_e, epsilon, time_step)
ps_ouss_asymptotic(freq, power_o, sigma, 
                   rho, lambda, 
                   power_e, epsilon, time_step)

Arguments

`freq`	Single number or numeric vector. The frequency for which to the power spectrum is to be calculated.
`power_o`	Single non-negative number. Power spectrum at zero-frequency generated by the underlying OU process, when sampled at the given `time_step`. Either `power_o` or `sigma` (but not both) must be provided.
`sigma`	Single number. Standard deviation of OU fluctuations around equilibrium. Either `power_o` or `sigma` (but not both) must be provided.
`rho`	Single number between 0 (exclusive) and 1 (inclusive). Correlation of the OU process between two subsequent time points. Either `rho` or `lambda` (but not both) must be provided.
`lambda`	Single non-negative number. Resilience (or relaxation rate) of the OU process. This is also the inverse correlation time of the OU process. Either `rho` or `lambda` (but not both) must be provided.
`power_e`	Single non-negative number. Asymptotic power spectrum at large frequencies due to the random measurement errors. Setting this to zero corresponds to the classical OU process. Either `power_e` or `epsilon` (but not both) must be provided.
`epsilon`	Single number. Standard deviation of Gaussian measurement error. Setting this to zero corresponds to the classical OU process. Either `power_e` or `epsilon` (but not both) must be provided.
`time_step`	Positive number. The time step of the time series that was (or will be) used for periodogram generation.

Details

The OUSS parameters power_o, lambda and power_e will typically be maximum-likelihood fitted values returned by evaluate.pm. time_step is also returned by evaluate.pm and is inferred from the analysed time series. More generally, power_o and power_e are proportional to the OUSS parameters sigma^2 and epsilon^2 (see generate_ouss), respectively, but the exact scaling depends on the normalization used for the periodogram.

Value

Returns a numeric vector of the same size as freq, containing the corresponding powers of the OUSS process.

Note

This function is the asymptotic version of ps_ouss in the limit where series_size becomes very large. If you want to compare the expected periodogram to the periodogram of a short time series use ps_ouss instead.

Author(s)

Stilianos Louca

References

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

Examples

# generate OUSS time series
times  = seq(0,20,0.25);
signal = generate_ouss(times, mu=0, sigma=1, lambda=1, epsilon=0.5);

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

# plot periodogram
plot(report$frequencies, report$periodogram, 
     type="l", ylab="power", xlab="frequency", main="periodogram & fitted OUSS power spectrum");

# plot OUSS power spectrum
lines(report$frequencies, 
      ps_ouss_asymptotic( freq=report$frequencies, 
                          power_o=report$power_o, 
                          lambda=report$lambda, 
                          power_e=report$power_e, 
                          time_step=report$time_step), 
      col="red");

# generate OUSS time series
times  = seq(0,20,0.25);
signal = generate_ouss(times, mu=0, sigma=1, lambda=1, epsilon=0.5);

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

# plot periodogram
plot(report$frequencies, report$periodogram, 
     type="l", ylab="power", xlab="frequency", main="periodogram & fitted OUSS power spectrum");

# plot OUSS power spectrum
lines(report$frequencies, 
      ps_ouss_asymptotic( freq=report$frequencies, 
                          power_o=report$power_o, 
                          lambda=report$lambda, 
                          power_e=report$power_e, 
                          time_step=report$time_step), 
      col="red");

Run example `evaluate.pm` analyses

Description

Generate random cyclic time series, analyse them using evaluate.pm and plot the results.

Usage

runExample(verbose=TRUE)
runExample(verbose=TRUE)

Arguments

verbose

Boolean, whether to print messages to the screen.

Value

No return value.

Author(s)

Stilianos Louca

References

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

Examples

  # this might take a few seconds
  runExample();
# this might take a few seconds
  runExample();

Statistical significance of local periodogram peaks

Description

Calculate statistical significance for a secondary periodogram peak (i.e. a non-global periodogram maximum), based on the null hypothesis of an OUSS process.

Usage

significanceOfLocalPeak(power_o, lambda, power_e, 
                        time_step, series_size,
                        Nfreq, peakFreq, peakPower)
significanceOfLocalPeak(power_o, lambda, power_e, 
                        time_step, series_size,
                        Nfreq, peakFreq, peakPower)

Arguments

`power_o`	Positive number. Power at zero-frequency stemming from the underlying OU process.
`lambda`	Positive number. Resilience (or relaxation rate) of the OU process, in inverse time units. This is also the inverse correlation time of the OU process.
`power_e`	Non-negative number. Asymptotic power at large frequencies due to random measurement errors. Setting this to zero corresponds to the classical OU process.
`time_step`	Positive number. The time step of the time series that was used to calculate the periodogram.
`series_size`	Positive integer. The size of the time series for which the periodogram peak was calculated.
`Nfreq`	The number of frequencies from which the local periodogram peak was picked. Typically equal to the number of frequencies in the periodogram.
`peakFreq`	Single number. The frequency of the focal peak.
`peakPower`	Single number. The periodogram power calculated for the focal peak.

Details

The OUSS parameters power_o, lambda and power_e will typically be maximum-likelihood fitted values returned by evaluate.pm. The time_step is also returned by evaluate.pm and is inferred from the analysed time series. The examined periodogram peak (as defined by peakFreq) will typically be a secondary peak of interest, masked by other stronger peaks or a low-frequency maximum. The significance of such a peak is not defined by standard tests.

Value

The returned P-value (referred to as “local P-value”) is the probability that an OUSS process with the specified parameters would generate a periodogram with a power-to-expectation ratio greater than peakPower/E, where E is the power spectrum of the OUSS process at frequency peakFreq. Hence, the significance is a measure for how much the peak power deviates from its expectation. The calculated value is an approximation. It becomes exact for long regular time series.

Note

This statistical significance is not equivalent to the one calculated by evaluate.pm for the global periodogram maximum. If the investigated periodogram peak is a global maximum, then the P-value returned by evaluate.pm should be preferred, as it also takes into account the absolute magnitude of the peak.

Author(s)

Stilianos Louca

References

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

Examples

  # In this example we generate a random cyclic time series, where the peak is (most likely)
  # masked by a strong low-frequency maximum.
  # We will use significanceOfLocalPeak() to evaluate its significance
  # based on its deviation from the expected power.
  
  # generate cyclic time series by adding a periodic signal to an OUSS process
  period    = 1;
  times     = seq(0,20,0.2);
  signal    = 0.5 * cos(2*pi*times/period) +
              generate_ouss(times, mu=0, sigma=1, lambda=1, epsilon=0.5);
  
  # calculate periodogram and fit OUSS model
  report    = evaluate.pm(times=times, signal=signal);
  print(report)
	
  # find which periodogram mode approximately corresponds to the frequency we are interested in
  cycleMode = which(report$frequencies>=0.99/period)[1]; 
	
  # calculate P-value for local 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]);

  # plot time series
  old.par <- par(mfrow=c(1, 2));
  plot(ts(times), ts(signal), 
       xy.label=FALSE, type="l", 
       ylab="signal", xlab="time", main="Time series (cyclic)", 
       cex=0.8, cex.main=0.9);

  # plot periodogram
  title = sprintf("Periodogram OUSS analysis\nfocusing on local peak at freq=%.3g\nPlocal=%.2g",
                  report$frequencies[cycleMode],Pvalue);
  plot(ts(report$frequencies), ts(report$periodogram), 
       xy.label=FALSE, type="l", 
       ylab="power", xlab="frequency", main=title, 
       col="black", cex=0.8, cex.main=0.9);
	
  # plot fitted OUSS power spectrum
  lines(report$frequencies, report$fittedPS, col="red");
	par(old.par)

# In this example we generate a random cyclic time series, where the peak is (most likely)
  # masked by a strong low-frequency maximum.
  # We will use significanceOfLocalPeak() to evaluate its significance
  # based on its deviation from the expected power.
  
  # generate cyclic time series by adding a periodic signal to an OUSS process
  period    = 1;
  times     = seq(0,20,0.2);
  signal    = 0.5 * cos(2*pi*times/period) +
              generate_ouss(times, mu=0, sigma=1, lambda=1, epsilon=0.5);
  
  # calculate periodogram and fit OUSS model
  report    = evaluate.pm(times=times, signal=signal);
  print(report)
	
  # find which periodogram mode approximately corresponds to the frequency we are interested in
  cycleMode = which(report$frequencies>=0.99/period)[1]; 
	
  # calculate P-value for local 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]);

  # plot time series
  old.par <- par(mfrow=c(1, 2));
  plot(ts(times), ts(signal), 
       xy.label=FALSE, type="l", 
       ylab="signal", xlab="time", main="Time series (cyclic)", 
       cex=0.8, cex.main=0.9);

  # plot periodogram
  title = sprintf("Periodogram OUSS analysis\nfocusing on local peak at freq=%.3g\nPlocal=%.2g",
                  report$frequencies[cycleMode],Pvalue);
  plot(ts(report$frequencies), ts(report$periodogram), 
       xy.label=FALSE, type="l", 
       ylab="power", xlab="frequency", main=title, 
       col="black", cex=0.8, cex.main=0.9);
	
  # plot fitted OUSS power spectrum
  lines(report$frequencies, report$fittedPS, col="red");
	par(old.par)

Package 'peacots'

Help Index

Periodogram peaks in correlated time series

Description

Details

Author(s)

References

Examples

Evaluate the periodogram maximum of a time series

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Statistical significance of periodogram peaks (classical)

Description

Usage

Arguments

Value

Author(s)

References

See Also

Examples

Generate random time series of the OUSS process

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Plot results of a evaluate.pm analysis

Description

Usage

Arguments

Value

Author(s)

References

See Also

Examples

Expected periodogram of the OUSS process

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Power spectrum of the OUSS process

Description

Usage

Arguments

Details

Value

Note

Author(s)

References

See Also

Examples

Run example evaluate.pm analyses

Description

Usage

Arguments

Value

Author(s)

References

See Also

Examples

Statistical significance of local periodogram peaks

Description

Usage

Arguments

Plot results of a `evaluate.pm` analysis

Run example `evaluate.pm` analyses