Description Usage Arguments Details Value Note Author(s) References See Also Examples
Calculate statistical significance for a secondary periodogram peak (i.e. a non-global periodogram maximum), based on the null hypothesis of an OUSS process.
1 2 3 | significanceOfLocalPeak(power_o, lambda, power_e,
time_step, series_size,
Nfreq, peakFreq, peakPower)
|
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. |
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.
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.
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.
Stilianos Louca
Louca, S., Doebeli, M. (2015) Detecting cyclicity in ecological time series, Ecology 96: 1724–1732
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 | # 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.