The standard approach to the application of spectral analysis in the search for orbitally forced cycles in stratigraphic data series involves comparing the power spectrum of the data with that of a statistical model representing (random) background noise. The discrimination of a (candidate) periodic cycle from random noise requires estimation of the statistical confidence that power at a given frequency in the spectrum is unlikely to be part of the noise. Andrews et al. (2016) made spectral analyses of natural gamma-ray (GR) data from onshore and offshore sections through the lacustrine Caithness Flags (Middle Devonian, Scotland), the offshore section comprising 807 m of strata in the UK Continental Shelf (UKCS) well 11/25-2 ST1. They generated periodograms using Thomson's multi-taper method (MTM), followed by the Matlab function Redconf (Husson et al. 2014) to model the noise in the data and to derive their statistical confidence levels. As is the convention in cyclostratigraphic analysis, Redconf uses a first-order autocorrelation (AR(1)) process to model the random noise inherent in the data. The resulting periodograms for the total dataset and for subgroup-level subsets (fig. 9a–e of Andrews et al. 2016) appear to show the presence of substantial numbers of cycles significant at the statistical 90% confidence level (CL).
This standard approach to the estimation of confidence levels in cyclostratigraphy was criticized by Vaughan et al. (2011). They found two sources of high false positive rates of cycle identification: (1) failure to correct the probabilistic CLs for the inevitably multiple nature of the significance tests; (2) universal reliance on AR(1) for modelling the noise. We show here that almost all of the cycles identified by Andrews et al. result from incorrect estimates of CLs.
Figure 1 is our bilogarithmic plot of the complete MTM periodogram of the full (5291-point) 11/25-2 ST1 Caithness Group GR dataset, detrended as in the original study, and with confidence limits calculated using Husson's Redconf. Compared with Andrews et al. (2016, fig. 9a), our plot shows the full spectrum and, by including the 50% CL (the median of the noise model), better illustrates the fit of the model to the data. Figure 1 shows that the Redconf AR(1) model is an imperfect fit to the spectrum of the data, a fault that Vaughan et al. (2011) showed is an important source of false positive cycle identifications in cyclostratigraphic analyses. The conspicuous power minimum centred on a wavelength of 0.46 m is clearly an artefact resulting from filtering, probably at an early stage in processing of the log. Its distorting effect on the noise model is a partial explanation for the latter's poor fit to the data, notably in the 1–10 m wavelength range, where Andrews et al. (2016, fig. 9a) found multiple cycles, although only two were selected for interpretation.
MTM periodogram of 11/25-2 ST1 GR data (all 5291 points) after subtraction of a third-order polynomial trend. The AR(1) noise model and confidence limits were estimated using the method of Husson et al. (2014). Compare figure 9a of Andrews et al. (2016), which shows only part of this spectrum.
In Figure 2, we have used the procedure adopted by Andrews et al. on a version of the 11/25-2 data from which two out of every three data-points have been removed. This reduction of the data eliminates the shortest wavelengths from the power spectrum, including the artificial minimum at 0.46 m; it also decreases the correlation between successive GR readings (see below). Other than raising the minimum wavelength that can be detected, it has no effect on the data's power spectrum. It improves the fit between the spectrum of the data and the median (the 50% CL) of the Redconf-calculated AR(1) noise model, but there are still multiple peaks exceeding the 90% CL, highlighting a second and more important source of false positive (FP) ‘cycle’ identifications.
MTM periodogram of 11/25-2 GR data (1764 points) after removing two of every three data-points and subtracting a third-order polynomial trend. Confidence limits at 50%, 90% and 99.99% (the last corresponding to a global confidence limit of 90%) were calculated as in Figure 1.
The numerous significant peaks in Figure 1 arise from incorrect estimation of the CLs. The CLs are statistical thresholds; spectral power exceeding the threshold at any single frequency is unlikely (at CL%) to belong to the spectrum of the chosen noise model. For a dataset of N points, the periodogram comprises estimates of spectral power at Nf frequencies, where Nf = N/2. Confidence limits as estimated by Redconf (and by most other implementations of the standard cyclostratigraphic method) are relevant to any significance test applied at a single frequency. Such tests ask ‘Is the spectral power at frequency F statistically distinct from the modelled noise spectrum at that frequency?’ In cyclostratigraphy, the time–depth relationship of the data is not known a priori, so it is unknown which spectral frequencies should be tested as candidate Milankovitch periods. Statistical significance must therefore be tested at every periodogram point: we ask instead ‘Are there any frequencies at which the data's spectral power is statistically unlikely to be part of the modelled noise spectrum?’. This requires a correction to the CLs in proportion to the number of frequencies in the periodogram; otherwise spectral power will appear, incorrectly, to exceed the local CL at a large number of frequencies, as it does in Figures 1 and 2.
Any CL entails a complementary rate of false positive cycle identifications: CL + FP = 1 (or CL% + FP% = 100%). The necessary correction for multiple tests can be achieved by dividing the desired global (whole-periodogram) FP rate by the number of periodogram frequencies (for details see Vaughan et al. 2011). Andrews et al. used a CL of 90% (0.9), which is a false positive rate of 10% (0.1). For the reduced 11/25-2 data (1762 points; 881 frequencies), the corrected single-test CL is therefore (1 − 0.1/881), or c. 99.99%. The periodogram can now be scanned for points that exceed this single-test CL (the dashed line in Fig. 2). No peaks in Figure 2 exceed this CL, meaning that, even at this relatively low (90%) level of statistical significance, there are no cycles present in the data. (The application of this ‘Bonferroni correction’ requires that the data-points, and hence the periodogram frequencies, are mutually independent. GR logging typically imposes a significant degree of correlation between adjacent data-points, but the omission of two out of every three GR values should largely have eliminated this.) We have applied the same approach to the four individually detrended subsets of the (reduced) Caithness Flags GR data (compare Andrews et al. 2016, fig. 9b–e). The Mey Subgroup alone yields a single periodogram peak that reaches the global 90% CL, which (for this dataset of 755 points) corresponds to a single-test CL of 99.97% (Fig. 3). This peak is at a frequency of 0.38 cycles m−1, a wavelength of 2.63 m, which corresponds to the shortest of the three ‘significant’ wavelengths found by Andrews et al. (their fig. 9b). The other three Caithness subgroups yielded no spectral peaks at this level of significance. The time calibration of the Caithness Group by Andrews et al. (2016) is of correspondingly doubtful value.
MTM periodogram for Mey Formation in 11/25-2 (reduced GR data, 755 data-points, 378 periodogram points) after subtracting a third-order polynomial trend. Confidence limits at 50%, local 90% and 99.97% (=global 90% CL).
Concerns have been expressed by, for example, Hilgen et al. (2015) that too stringent a multiple-test correction risks ‘false negative’ results: real cycles in the data are missed because the CLs are raised too high. Kemp (2016), agreeing that a correction is unavoidable, has started to explore the balance of risks between false positive and false negative results with experiments on data that conceal a cycle period in random red noise. For cyclostratigraphic investigations where there can be no certainty of the existence of any ‘real’ cycle period(s), we continue to urge caution.
Acknowledgements
We thank S. Andrews and D. Cornwell for providing their data and introducing us to Husson's Matlab code, and A. Johnson for correspondence about the effects of filtering on well logs. We also thank Prof. D. Waltham for his constructive review.
Scientific editing by Stuart Jones
- © 2018 The Author(s)
Discussion on ‘A 2.3 million year lacustrine record of orbital forcing from the Devonian of northern Scotland’, Journal of the Geological Society, London 173, 474 – 488
