lowspec {astrochron}

Robust Locally-Weighted Regression Spectral Background Estimation

Description

LOWSPEC: Robust Locally-Weighted Regression Spectral Background Estimation (Meyers, 2012)

Usage

lowspec(dat,decimate=NULL,tbw=3,padfac=5,detrend=F,siglevel=0.9,setrho,
         lowspan,b_tun,output=0,CLpwr=T,xmin,xmax,pl=1,sigID=T,genplot=T,
         verbose=T)

Arguments

dat

Stratigraphic series for LOWSPEC. First column should be location (e.g., depth), second column should be data value.

decimate

Decimate statigraphic series to have this sampling interval (via piecewise linear interpolation). By default, no decimation is performed.

tbw

MTM time-bandwidth product (2 or 3 permitted)

padfac

Pad with zeros to (padfac*npts) points, where npts is the original number of data points.

detrend

Remove linear trend from data series? This detrending is performed following AR1 prewhitening. (T or F)

siglevel

Significance level for peak identification. (0-1)

setrho

Define AR1 coefficient for pre-whitening (otherwise calculated). If set to 0, no pre-whitening is applied.

lowspan

Span for LOWESS smoothing of prewhitened signal, usually fixed to 1. If using value <1, the method is overly conservative with a reduced false positive rate.

b_tun

Robustness weight parameter for LOWSPEC. By default, this will be estimated internally.

output

What should be returned as a data frame? (0=nothing; 1=pre-whitened spectrum + harmonic F-test CL + LOWSPEC background + LOWSPEC CL + 90%-99% LOWSPEC power levels; 2=sig peaks)

CLpwr

Plot LOWSPEC noise confidence levels on power spectrum? (T or F)

xmin

Smallest frequency for plotting.

xmax

Largest frequency for plotting.

pl

Power spectrum plotting: (1) linear frequency-log spectral power, (2) linear frequency-linear spectral power (3) log frequency-log spectral power, (4) log frequency-linear spectral power

sigID

Identify significant frequencies on power and probabilty plots? (T or F)

genplot

Generate summary plots? (T or F)

verbose

Verbose output? (T or F)

Details

LOWSPEC is a 'robust' method for spectral background estimation, designed for the identification of potential astronomical signals that are imbedded in red noise (Meyers, 2012). The complete algoritm implemented here is as follows: (1) initial pre-whitening with AR1 filter (default) or other filter as appropriate (e.g., see function prewhiteAR), (2) power spectral estimation via the multitaper method (Thomson, 1982), (3) robust locally weighted estimation of the spectral background using the LOWESS-based (Cleveland, 1979) procedure of Ruckstuhl et al. (2001), (4) assignment of confidence levels using a Chi-square distribution.

NOTE: If you choose to pre-whiten before running LOWSPEC (rather than using the default AR1 pre-whitening), specify setrho=0.

Candidiate astronomical cycles are subsequently identified via isolation of those frequencies that achieve the required (e.g., 90 percent) LOWSPEC confidence level and MTM harmonic F-test confidence level. Allowance is made for the smoothing inherent in the MTM power spectral estimate as compared to the MTM harmonic spectrum. That is, an F-test peak is reported if it achieves the required MTM harmonic confidence level, while also achieving the required LOWSPEC confidence level within +/- half the power spectrum bandwidth resolution. One additional criterion is included to further reduce the false positive rate, a requirement that significant F-tests must occur on a local power spectrum high, which is parameterized as occurring above the local LOWSPEC background estimate. See Meyers (2012) for futher information on the algorithm.

In this implementation, the 'robustness criterion' ('b' in EQ. 6 of Ruckstuhl et al., 2001) has been optimized for 2 and 3 pi DPSS, using a 'span' of 1. By default the robustness criterion will be estimated. Both 'b' and the 'span' can be expliclty set using parameters 'b_tun' and 'lowspan'. Note that it is permissible to decrease 'lowspan' from its default value, but this will result in an overly conservative false positive rate. However, it may be necessary to reduce 'lowspan' to provide an approporiate background fit for some stratigraphic data. Another option is to decimate the data series prior to spectral estimation.

Value

If option 1 is selected, a data frame containing the following is returned: Frequency, Prewhitened power, harmonic F-test CL, LOWSPEC CL, LOWSPEC background, 90%-99% LOWSPEC power levels. NOTE: as of version 0.8, the order of the columns in the output data frame has been changed, for consistency with functions mtm, mtmML96, and mtmPL.

If option 2 is selected, the 'significant' frequencies are returned (as described above).

References

W.S. Cleveland, 1979, Locally weighted regression and smoothing scatterplots: Journal of the American Statistical Association, v. 74, p. 829-836.

S.R. Meyers, 2012, Seeing Red in Cyclic Stratigraphy: Spectral Noise Estimation for Astrochronology: Paleoceanography, 27, PA3228, doi:10.1029/2012PA002307.

A.F. Ruckstuhl, M.P Jacobson, R.W. Field, and J.A. Dodd, 2001, Baseline subtraction using robust local regression estimation: Journal of Quantitative Spectroscopy & Radiative Transfer, v. 68, p. 179-193.

D.J. Thomson, 1982, Spectrum estimation and harmonic analysis: IEEE Proceedings, v. 70, p. 1055-1096.

See Also

eha, mtm, mtmAR, mtmML96, periodogram, rfbaseline, and spec.mtm

Examples

# generate example series with periods of 400 ka, 100 ka, 40 ka and 20 ka
ex = cycles(freqs=c(1/400,1/100,1/40,1/20),start=1,end=1000,dt=5)

# add AR1 noise
noise = ar1(npts=200,dt=5,sd=.5)
ex[2] = ex[2] + noise[2]

# LOWSPEC analysis
pl(1, title="lowspec")
lowspec(ex)

# compare to MTM spectral analysis, with conventional AR1 noise test
pl(1,title="mtm")
mtm(ex,ar1=TRUE)

# compare to ML96 analysis
pl(1, title="mtmML96")
mtmML96(ex)

# compare to amplitudes from eha
pl(1,title="eha")
eha(ex,tbw=3,win=1000,pad=1000)