| Title: | Lomb-Scargle Periodogram |
|---|---|
| Description: | Computes the Lomb-Scargle Periodogram and actogram for evenly or unevenly sampled time series. Includes a randomization procedure to obtain exact p-values. Partially based on C original by Press et al. (Numerical Recipes) and the Python module Astropy. For more information see Ruf, T. (1999). The Lomb-Scargle periodogram in biological rhythm research: analysis of incomplete and unequally spaced time-series. Biological Rhythm Research, 30(2), 178-201. |
| Authors: | Thomas Ruf [aut, cre] |
| Maintainer: | Thomas Ruf <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 2.5.0 |
| Built: | 2024-11-25 06:26:16 UTC |
| Source: | https://github.com/thomaspruf/lomb |
The Lomb-Scargle periodogram is the most widely used method to detect even weak periodic components in unequally sampled time series. It can also be used for equally sampled time series.The oacka alao plots actograms and converts them to periodograms.
| Package: | lomb |
| Type: | Package |
| Version: | 2.1.0 |
| Date: | 2022-02-22 |
| License: | GPL-3 |
Function lsp computes the Lomb-Scargle periodogram for unevenly sampled times series (e.g., series with missing data). P-values for the highest peak in the periodogram are computed from the exponential distribution. Alternatively, function randlsp computes a p-value for the largest peak in the periodogram by repeatedly randomising the time-series sequence. Both functions allow setting the range of frequencies to be inspected, as well as the stepsize (oversampling factor) used for frequency scanning. Function actogram plots actograms and function makedf prepares them for lsp.
Thomas Ruf
Department of Interdisciplinary Life Sciences,
University of Veterinary Medicine, Vienna, Austria
Maintainer: Thomas Ruf [email protected]
Ruf, T. (1999) The Lomb-Scargle Periodogram in Biological Rhythm Research: Analysis of Incomplete and Unequally Spaced Time-Series. Biological Rhythm Research 30: 178–201
data(lynx) lsp(lynx)data(lynx) lsp(lynx)
plots an actogram for a time series with irregular (or regular) sampling intervals.
actogram(date, response,from ,to, scalefac, subtract, dble, dig, border, fill, grad, lwd, photo,latitude, longitude, zone, twilight)actogram(date, response,from ,to, scalefac, subtract, dble, dig, border, fill, grad, lwd, photo,latitude, longitude, zone, twilight)
date |
data datetime, format as in as.Date |
response |
which variable to plot |
from |
date for start of subsection |
to |
date for end of subsection |
scalefac |
one day is 1.0 wide. Use >1 if plots are exaggerated |
subtract |
value to subtract from response. |
dble |
logical Double plot? Ie day1 day2, day2 day3, day3 day4 ... |
dig |
logical Digitize to 0 1 based on mean? |
border |
border_colour of rectangle edges |
fill |
colour of rectangle fills, relevant only if timepoint separation is large |
grad |
logical Plot gradient? |
lwd |
line width of rectangles |
photo |
logical Plot photoperiod? Photo is true when location is entered. |
latitude |
animal location |
longitude |
animal location |
zone |
time zone of location. Default:0 = Greenwich |
twilight |
"rise/set", "civil" or "nautic" |
This function plots actograms of both unevenly and evenly sampled data. It only requires data in standard R format, no special rhythms format is necessary.
## Not run: data(caradat) actogram(caradat$Date, caradat$Activity, dble=TRUE, photo=FALSE, dig=TRUE, fill="blue") data(deerdat) actogram(deerdat$Zeit, deerdat$Akt, grad=TRUE,from="2010-10-01 00:00:00",to="2011-03-31 00:00:00", latitude=47.1415,longitude=9.5215,zone=1,twilight="nautic") ## End(Not run)## Not run: data(caradat) actogram(caradat$Date, caradat$Activity, dble=TRUE, photo=FALSE, dig=TRUE, fill="blue") data(deerdat) actogram(deerdat$Zeit, deerdat$Akt, grad=TRUE,from="2010-10-01 00:00:00",to="2011-03-31 00:00:00", latitude=47.1415,longitude=9.5215,zone=1,twilight="nautic") ## End(Not run)
Locomotor activiy of a blind beetle.
data("caradat")data("caradat")
A data frame with 2014 observations on the following 2 variables.
a numeric vector with date and time of day
a numeric vector of locomotor activity
activity under DD in a cave observed at UNEVEN intervals.
The data were kindly provided by F. Weber, Münster, Germany. All the experiments were carried out between 1973 and 1980. The dates (but not their order and times) in the file are fictitious.
Locomotor activiy and other variables of a red deer free-living in the alps.
data("deerdat")data("deerdat")
A data frame with 293826 observations on the following 9 variables.
idanimal ID
Tiernummeragain
Halsbandnummercollar number
Zeittime
RepeaterTempgevice temperature
Bodytempbody temperature
Kopfwechselnumber of head down movements
KopfuntenzeitSektime head down
Aktactivity
data(deerdat) ## maybe str(deerdat) ; plot(deerdat) ...data(deerdat) ## maybe str(deerdat) ; plot(deerdat) ...
Retrieves and displays the npeaks largest peaks in the periodogram-
getpeaks(object,npeaks,plotit)getpeaks(object,npeaks,plotit)
object |
object must be of class "lsp" |
npeaks |
number of peaks to get |
plotit |
if TRUE show plot |
Returns a list with
data |
A dataframe with times an heights of peaks |
plot |
An annotated periodogram |
Thomas Ruf [email protected]
per=lsp(lynx,ofac=5) getpeaks(per,6) # obtain the 6 largest peaksper=lsp(lynx,ofac=5) getpeaks(per,6) # obtain the 6 largest peaks
From astropy.timeseries
ggamma(N)ggamma(N)
N |
A positive number |
sqrt(2 / N) * exp(lgamma(N / 2) - lgamma((N - 1) / 2))
Thomas Ruf [email protected]
VanderPlas, J. & Ivezic, Z. (2015) Periodograms for Multiband Astronomical Time Series.The Astrophysical Journal 812.1:18
ggamma(3)ggamma(3)
Telemetric measurements of rumen temperature in a free-living alpine ibex (Capra ibex) measured at unequal time intervals.
data(ibex)data(ibex)
A data frame with 1201 observations on 3 variables.
a character variable giving date and time of measurements.
a numerical variable giving hours elapsed since the first measurement.
a numerical variable giving rumen (stomach) temperature in degrees Celsius.
A subset of data from Signer, C., Ruf, T., Arnold, W. (2011) Functional Ecology 25: 537-547.
data(ibex) datetime <- as.POSIXlt(ibex$date) plot(datetime,ibex$temp,pch=19,cex=0.3)data(ibex) datetime <- as.POSIXlt(ibex$date) plot(datetime,ibex$temp,pch=19,cex=0.3)
activity and body temperature of a domestic dog in summer
data("layla")data("layla")
A data frame with 10120 observations on the following 4 variables.
a vector with animal ID
a vector with date and time of day
a numeric vector of body temperature
a numeric vector of kocomtor activity
data(layla) ## maybe str(layla) ; plot(layla) ...data(layla) ## maybe str(layla) ; plot(layla) ...
utility function to determine deviation from p-value
levopt(x, alpha, fmax, tm)levopt(x, alpha, fmax, tm)
x |
vector with start values |
alpha |
desired significance level |
fmax |
the highest frequency inspected |
tm |
a vector with measurement timepoints |
(log(prob)-log(alpha))^2
Thomas Ruf [email protected]
Computes the Lomb-Scargle periodogram for a time series with irregular (or regular) sampling intervals. Allows selecting a frequency range to be inspected, as well as the spacing of frequencies scanned.
lsp(x, times = NULL, from = NULL, to = NULL, type = c("frequency", "period"), ofac = 1, alpha = 0.01, normalize=c("standard","press"), plot = TRUE, ...)lsp(x, times = NULL, from = NULL, to = NULL, type = c("frequency", "period"), ofac = 1, alpha = 0.01, normalize=c("standard","press"), plot = TRUE, ...)
x |
The data to be analysed. x can be either a two-column numerical dataframe or matrix, with sampling times in column 1 and measurements in column 2, a single numerical vector containing measurements, or a single vector |
times |
If x is a single vector, times can be provided as a numerical vector of equal length containing sampling times. If x is a vector and times is NULL, the data are assumed to be equally sampled and times is set to 1:length(x). |
from |
The starting frequency (or period, depending on type) to begin scanning for periodic components. |
to |
The highest frequency (or period, depending on type) to scan. |
type |
Either “frequency” (the default) or “period”. Determines the type of the periodogram x-axis. |
ofac |
The oversampling factor. Must be an integer>=1. Larger values of ofac lead to finer scanning of frequencies but may be time-consuming for large datasets and/or large frequency ranges (from...to). |
alpha |
The significance level. The periodogram plot shows a horizontal dashed line. Periodogram peaks exceeding this line can be considered significant at alpha. Defaults to 0.01. Only used if plot=TRUE. |
normalize |
The type of normalization used, either “standard” or “press”. If normalization is standard (the default) the periodogram is confined to the interval 0-1, and the statistical significance of the largest peak in the periodogram is computed according to Baluev (2008).if normalization is set to “press” the periodogram will be normalized using the factor 1/(2 * var(y)) and the p-value for the significance of the largest peak in the periodogram is computed from the exponential distribution, as outlined in Press et al. (1994), see below |
plot |
Logical. If plot=TRUE the periodogram is plotted. |
... |
Further graphical parameters affecting the periodogram plot. |
For a more robust - but potentially time-consuming estimation of p-values (when n is large) see randlsp.
Significance levels in both lsp and randlsp increase with the number of frequencies inspected. Therefore, if the frequency-range of interest can be narrowed down a priori, use arguments “from” and “to” to do so.
A named list with the following components:
normalize |
The type of normalization used. |
scanned |
A vector containing the frequencies/periods scanned. |
power |
A vector containing the normalised power corresponding to scanned frequencies/periods. |
data |
Names of the data vectors analysed. |
n |
The length of the data vector. |
type |
The periodogram type used, either "frequency" or "period". |
ofac |
The oversampling factor used. |
n.out |
The length of the output (powers). This can be >n if ofac >1. |
alpha |
The false alarm probability used. |
sig.level |
Powers > sig.level can be considered significant peaks at p=alpha. |
peak |
The maximum power in the frequency/period interval inspected. |
peak.at |
The frequency/period at which the maximum peak occurred. |
p.value |
The probability that the maximum peak occurred by chance. |
For a description of the properties of the Lomb-Scargle Periodogram, its computation and comparison with other methods see Ruf, T. (1999). Function lsp uses the algorithm given by Press et al (1994). The Lomb-Scargle Periodogram was originally proposed by Lomb N.R. (1976) and further extended by Scargle J.D. (1982). An improved method for assessing the statistical significance of candidate periodicities by Baluev (2008), based on extreme value theory, is also implemented. This implementation uses code modified from the astropy.timeseries Python package (VanderPlas et al. 2012, 2015).
Thomas Ruf [email protected] based on code by Press et al (1994).
Baluev, R. V. (2008). Assessing the statistical significance of periodogram peaks. Monthly Notices of the Royal Astronomical Society, 385(3), 1279-1285.
Lomb N.R. (1976) Least-squares frequency analysis of unequally spaced data. Astrophysics and Space Science 39:447–462
Press W.H., Teukolsky S.A., Vetterling S.T., Flannery, B.P. (1994) Numerical recipes in C: the art of scientific computing.2nd edition. Cambridge University Press, Cambridge, 994pp.
Ruf, T. (1999) The Lomb-Scargle Periodogram in Biological Rhythm Research: Analysis of Incomplete and Unequally Spaced Time-Series. Biological Rhythm Research 30: 178–201.
Scargle J.D. (1982) Studies in astronomical time series. II. Statistical aspects of spectral analysis of unevenly spaced data. The Astrophysical Journal 302: 757–763.
VanderPlas, J., Connolly, A. Ivezic, Z. & Gray, A. (2012) Introduction to astroML: Machine learning for astrophysics. Proceedings of the Conference on Intelligent Data Understanding
VanderPlas, J. & Ivezic, Z. (2015) Periodograms for Multiband Astronomical Time Series.The Astrophysical Jounal 812.1:18
# ibex contains an unevenly sampled time series data(ibex) lsp(ibex[,2:3],ofac=5) lsp(ibex$temp,times=ibex$hours,type='period',ofac=5) # lynx contains evenly sampled data lsp(lynx) lynx.spec <- lsp(lynx,type='period',from=2,to=20,ofac=5) summary(lynx.spec) # generate unevenly sampled data time=(runif(200,1,1000)) y=2*cos(time/6)+rnorm(200,0,4) lsp(y,times=time,ofac=10, to=0.3)# ibex contains an unevenly sampled time series data(ibex) lsp(ibex[,2:3],ofac=5) lsp(ibex$temp,times=ibex$hours,type='period',ofac=5) # lynx contains evenly sampled data lsp(lynx) lynx.spec <- lsp(lynx,type='period',from=2,to=20,ofac=5) summary(lynx.spec) # generate unevenly sampled data time=(runif(200,1,1000)) y=2*cos(time/6)+rnorm(200,0,4) lsp(y,times=time,ofac=10, to=0.3)
Converts an actogram to a periodogram
makedf(tvar, pvar)makedf(tvar, pvar)
tvar |
data datetime |
pvar |
which variable to plot |
a data.frame with two colums: time and variable (eg. activity)
## Not run: data(caradat) #unevenly sampled focus=actogram(caradat$Date, caradat$Activity, dble=TRUE, photo=FALSE, zone=1, from="1970-01-01 00:00:00",to="1970-01-14 00:00:00") df=makedf (focus$date, focus$plotvar) lsp(df, type="period",ofac=5,from=12,to=36) data(layla) #evenly sampled focus=actogram(layla$Date,layla$Activity,latitude=48.20, longitude=16.37, zone=2, dig=TRUE) df=makedf (focus$date, focus$plotvar) lsp(df, type="period",ofac=5,from=20,to=50) ## End(Not run)## Not run: data(caradat) #unevenly sampled focus=actogram(caradat$Date, caradat$Activity, dble=TRUE, photo=FALSE, zone=1, from="1970-01-01 00:00:00",to="1970-01-14 00:00:00") df=makedf (focus$date, focus$plotvar) lsp(df, type="period",ofac=5,from=12,to=36) data(layla) #evenly sampled focus=actogram(layla$Date,layla$Activity,latitude=48.20, longitude=16.37, zone=2, dig=TRUE) df=makedf (focus$date, focus$plotvar) lsp(df, type="period",ofac=5,from=20,to=50) ## End(Not run)
Computes the statistical significance of peaks (range 0-1) in the standardized perodogram. Typically not called by the user.
pbaluev(Z,fmax,tm)pbaluev(Z,fmax,tm)
Z |
the height of a periodogram peak |
fmax |
the highest frequency inspected |
tm |
a vector with measurement timepoints |
Based on results in extreme value theory, improved analytic estimations of false alarm probabilities are given.
Returns the significance of the largest peak in the periodogram.
Code based on astropy.timeseries
Thomas Ruf [email protected].
Baluev, R. V. (2008). Assessing the statistical significance of periodogram peaks. Monthly Notices of the Royal Astronomical Society, 385(3), 1279-1285.
pbaluev(0.19,2.0,1:100)pbaluev(0.19,2.0,1:100)
Shows a periodogram in browser window as line and dot plot. When moving the cursor close to dots times an peak-heights of the periodogram are shown.
pershow(object) # object of class "lsp"pershow(object) # object of class "lsp"
object |
an object of class "lsp" |
Thomas Ruf [email protected]
per=lsp(lynx,ofac=2) pershow(per) #In Rstudio go to the viewer pane. Move mouse to point of interest.per=lsp(lynx,ofac=2) pershow(per) #In Rstudio go to the viewer pane. Move mouse to point of interest.
computes sunrise & sunset for day of year and location
photoperiod(dayofyear, latitude, longitude, zone = 0, twilight = "civil")photoperiod(dayofyear, latitude, longitude, zone = 0, twilight = "civil")
dayofyear |
day of ywar 1-366 |
latitude |
e.g. 42.00 |
longitude |
e.g. 9.00 |
zone |
time zone e.g 1 (Vienna) |
twilight |
"rise/set", "civil" or "nautic" |
sunrise |
vector of sunrises |
set |
vector of sunsets |
http://lexikon.astronomie.info/zeitgleichung/
photoperiod (180, 42,9, zone=1, twilight="civil")photoperiod (180, 42,9, zone=1, twilight="civil")
Plots the normalised power as a function of frequency (or period, depending on type in function lsp).
## S3 method for class 'lsp' plot(x, main = "Lomb-Scargle Periodogram", xlabel = NULL, ylabel = "normalized power", level = TRUE, plot=TRUE, ...)## S3 method for class 'lsp' plot(x, main = "Lomb-Scargle Periodogram", xlabel = NULL, ylabel = "normalized power", level = TRUE, plot=TRUE, ...)
x |
Object of class lsp as returned from function lsp. |
main |
Character. Main title of the periodogram plot. Defaults to “Lomb-Sargle Periodogram”. |
xlabel |
Character. X-axis label of the periodogram plot. |
ylabel |
Character. Y-axis label of the periodogram plot. |
level |
Logical. If TRUE, the significance level is displayed as a dashed line. |
plot |
If TRUE, the periodogram is plotted. |
... |
Additional graphics parameters |
Usually, this function is only called by function lsp. It maybe called by the user for some control of the output. For better control, plot results from lsp ($scanned, $power) as desired.
Invisibly returns the object of class lsp it is called with.
Thomas Ruf [email protected]
data(ibex) ibex.spec <- lsp(ibex[,2:3],type='period', from=12,to=36,ofac=10, plot=FALSE) plot.lsp(ibex.spec, main="Tb in Capra ibex",xlabel="Period (h)",ylabel="Power",level=FALSE)data(ibex) ibex.spec <- lsp(ibex[,2:3],type='period', from=12,to=36,ofac=10, plot=FALSE) plot.lsp(ibex.spec, main="Tb in Capra ibex",xlabel="Period (h)",ylabel="Power",level=FALSE)
randlsp is used to obtain robust p-values for the significance of the largest peak in a Lomb-Scargle periodogram by randomisation. The data sequence is scrambled repeatedly and the probability of random peaks reaching or exceeding the peak in the original (unscrambled) periodogram is computed.
randlsp(repeats=1000,x, times = NULL, from = NULL, to = NULL, type = c("frequency", "period"), ofac = 1, alpha = 0.01, plot = TRUE, trace = TRUE, ...)randlsp(repeats=1000,x, times = NULL, from = NULL, to = NULL, type = c("frequency", "period"), ofac = 1, alpha = 0.01, plot = TRUE, trace = TRUE, ...)
repeats |
An integer determining the number of repeated randomisations. Large numbers (>=1000) are better but can make the procedure time-consuming. |
x |
The data to be analysed. x can be either a two-column numerical dataframe or matrix, with sampling times in columnn 1 and measurements in column 2, a single numerical vector containing measurements, or a single vector |
times |
If x is a single vector, times can be provided as a numerical vector of equal length containing sampling times. If x is a vector and times is NULL, the data are assumed to be equally sampled and times is set to 1:length(x). |
from |
The starting frequency (or period, depending on type) to begin scanning for periodic components. |
to |
The highest frequency (or period, depending on type) to scan. |
type |
Either “frequency” (the default) or “period”. Determines the type of the periodogram x-axis. |
ofac |
The oversampling factor. Must be an integer >=1. Larger values of ofac lead to finer scanning of frequencies but may be time-consuming for large datasets and/or large frequency ranges (from...to). |
alpha |
The significance level. The periodogram plot shows a horizontal dashed line. Periodogram peaks exceeding this line can be considered significant at alpha. Defaults to 0.01. Only used if plot=TRUE. |
plot |
Logical. If TRUE, two plots are displayed (i) The periodogram of the original (unscrambled) data (ii) A histogram of peaks occurring by chance during sequence randomisation. A vertical line is drawn at the height of the peak in a periodogram of the original data. |
trace |
Logical. If TRUE, information about the progress of the randomisation procedure is printed during the running of randlsp. |
... |
Additional graphical parameters affecting the histogram plot. |
Function randlsp preserves the actual measurement intervals, which may affect the periodogram (see Nemec & Nemec 1985, below). Hence, this is a conservative randomisation procedure.
P-values from both randlsp and lsp increase with the number of frequencies inspected. Therefore, if the frequency-range of interest can be narrowed down a priori, use arguments “from” and “to” to do so.
A named list with the following items:
scanned |
A vector containing the frequencies/periods scanned. |
power |
A vector containing the normalised power corresponding to scanned frequencies/periods. |
data |
Names of the data vectors analysed. |
n |
The length of the data vector. |
type |
The periodogram type used, either “frequency” or “period”. |
ofac |
The oversampling factor used. |
n.out |
The length of the output (powers). This can be >n if ofac >1. |
peak |
The maximum power in the frequency/period interval inspected. |
peak.at |
The frequency/period at which the maximum peak occurred. |
random.peaks |
A vector of peaks (with length=repeats) of maximum power values computed from randomised data. |
repeats |
The number of randomisations. |
p.value |
The probability that the peak in the original data occurred by chance, computed from randomising the data sequence. |
Thomas Ruf [email protected]
Nemec A.F.L, Nemec J.M. (1985) A test of significance for periods derived using phase-dispersion-miminimization techniques. The Astronomical Journal 90:2317–2320
data(lynx) set.seed(444) rand.times <- sample(1:length(lynx),30) # select a random vector of sampling times randlsp(repeats=1000,lynx[rand.times],times=rand.times)data(lynx) set.seed(444) rand.times <- sample(1:length(lynx),30) # select a random vector of sampling times randlsp(repeats=1000,lynx[rand.times],times=rand.times)
Summary method for class lsp.
## S3 method for class 'lsp' summary(object,...)## S3 method for class 'lsp' summary(object,...)
object |
an object of class lsp. |
... |
currently, no other arguments are required. |
summary.lsp returns a one column data.frame with results from function lsp. Row names and contents are as follows:
Time |
Name of the sampling time variable. |
Data |
Name of the measured variable. |
Type |
either “frequency” or “period”. |
Oversampling factor |
The degree of oversampling (>=1). |
From |
The lowest frequency (or period, depending on type) inspected. |
To |
The highest frequency (or period, depending on type) inspected. |
# frequencies |
The number of frequencies (or periods, depending on type) inspected. |
PNmax |
The peak normalised power in the periodogram. |
At frequency |
The frequency at which PNmax occurred. |
At period |
The period at which PNmax occurred. |
P-value (PNmax) |
The probability that PNmax occurred by chance, computed from the exponential distribution. |
Thomas Ruf [email protected]
data(lynx) summary(lsp(lynx))data(lynx) summary(lsp(lynx))
Summary method for class randlsp.
## S3 method for class 'randlsp' summary(object,...)## S3 method for class 'randlsp' summary(object,...)
object |
an object of class randlsp. |
... |
currently, no other arguments are required. |
summary.randlsp returns a one column data.frame with results from function randlsp. Row names and contents are as follows:
Time |
Name of the sampling time variable. |
Data |
Name of the measured variable. |
Type |
either “frequency” or “period”. |
Oversampling |
The degree of oversampling (>=1). |
From |
The lowest frequency (or period, depending on type) inspected. |
To |
The highest frequency (or period, depending on type) inspected. |
# frequencies |
The number of frequencies (or periods, depending on type) inspected. |
PNmax |
The peak normalised power in the periodogram. |
At frequency |
The frequency at which PNmax occurred. |
At period |
The period at which PNmax occurred. |
Repeats |
The number of randomisations. |
P-value (PNmax) |
The probability that PNmax occurred by chance, computed from randomising the data sequence. |
Thomas Ruf [email protected]
data(lynx) summary(randlsp(repeats=500,lynx))data(lynx) summary(randlsp(repeats=500,lynx))
Import lsp ggplot2 theme. It builds on theme_bw.
theme_lsp(bs=18)theme_lsp(bs=18)
bs |
basesize of font |
A theme element
plot(lsp(lynx))+theme_lsp(25)plot(lsp(lynx))+theme_lsp(25)