peacots (version 1.0)

significanceOfLocalPeak: 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, 
                        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.

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.

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.

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. 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.

References

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

See Also

evaluate.pm

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.75 * 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, 
                                      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
focusing on local peak at freq=%.3g
Plocal=%.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)