significanceOfLocalPeak {peacots} | R Documentation |
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)
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
See Also
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)