Names the assumption behind a spectral fit. A periodogram ordinate of
a stationary Gaussian series is exponential about the spectral
density, so with a log link the linear predictor IS the log spectrum
and nothing else in the fit changes. whittle() is that exponential
family, with the periodogram-specific refusals attached.
Arguments
- tapers
Number of independent periodograms averaged into each ordinate.
1(the default) is the raw periodogram and gives the exponential likelihood;k > 1gives Gamma with the shape fixed atk.- link
Link for the mean. The default
"log"is what makes the linear predictor the log spectrum, and is very nearly always what is wanted.
Details
tapers = k for k > 1 is the Gamma likelihood with the shape held
at k, which is the distribution of an average of k independent
ordinates: a Welch or Bartlett estimate from
frm_periodogram(segments = k), or a k-taper multitaper estimate.
The shape is data preparation, not a parameter, so it is fixed rather
than estimated - through the same mechanism as bf(y ~ x, shape = k),
and an explicit shape in bf() still wins.
Post-processing speaks about the periodogram
fitted() returns the fitted spectral density at each row's
frequency, predict() the same for new frequencies, and
residuals() compares an ordinate with that density. None of them
is about the series, and none of them can be: the phases are not in
a periodogram.
simulate() is refused by name rather than documented, and so
are the four entry points built on it: pp_check(),
frm_simulate(), frm_bootstrap(), and
conditional_effects(method = "predict").
Draws from this family are ordinates, a vector of them under the
response's name reads as a simulated series, and that is a mistake a
user makes once and never notices. Each refusal names the two things
to do instead: draw ordinates in the open, with
rexp(nobs(fit), 1 / fitted(fit)) (or rgamma(nobs(fit), shape = k,
scale = fitted(fit) / k) at tapers = k), which is also the
one-line replacement for a DHARMa-style check; or draw a series
whose spectrum is the fitted one with frm_series_draw().
What the family refuses, and what it cannot see
Two things are checked, and both look only at the response.
A response that is not positive is refused outright: a periodogram ordinate is a squared modulus, so a negative one means log power or power in dB, which this likelihood is not about.
A response too SMOOTH to have the declared shape is also refused.
var(diff(log(y), lag = 3)) is 2 trigamma(tapers) for ordinates of
that shape whatever the spectrum is, and only the spectrum's own
variation across the step adds to it, so a value far below that is
evidence the ordinates were averaged more than tapers says, or
that they are spectral leakage rather than signal. The threshold is
a calibrated fraction of the model value that shrinks as the number
of ordinates falls.
The step is three ordinates rather than one because of the taper. A Hann window's transform is three bins wide, so it correlates NEIGHBORING ordinates (about 0.31 on the log scale) and leaves ordinates three apart uncorrelated (-0.003). At lag one a legitimate Hann-tapered response was refused about 2% of the time; at lag three the measured rate is 0 in 2000 replicates at each of nine cells, and the same holds for Hamming and Blackman windows, which this package does not apply and which were refused 1.0% and 18.1% of the time at lag one. The wider step costs leakage detection at the shortest usable length, 1.3 points at exponent 3 and 2.3 at exponent 2.5 with 127 ordinates, and under one point at 255 ordinates and above (2000 replicates a cell).
What it cannot see. It needs the rows in frequency order, which
is what frm_periodogram() returns: a periodogram sorted by its own
power is refused every time, wrongly. It does not run at all below
24 ordinates. It separates tapers = 2 from a raw periodogram only
above about 200 ordinates, because Gamma(2) and Gamma(1) are close.
It reads an estimate smoothed ACROSS FREQUENCY (a Daniell window, a
multitaper) poorly at both ends. Declared raw, such an estimate is
caught every time above about 100 ordinates and often missed below
it: a three-bin Daniell smooth is caught 18% of the time at 32
ordinates and 58% at 64. Declared honestly, by giving tapers the
equivalent degrees of freedom, it passes up to about five bins of
smoothing and is then refused anyway, 37% to 89% of the time at
seven bins and 92% or more at eleven, because the check assumes
ordinates independent across frequency and a wide smooth is not. Fit
a
frequency-smoothed estimate with care, or fit the raw or
segment-averaged periodogram it came from.
And it is a backstop, not a test to rely on: at 127 ordinates (a
one-second epoch at 256 Hz) an untapered exponent-3 power law is
caught about 92% of the time and an exponent-2 one much less often.
Choose the taper from the shape of the spectrum, not from whether an
error appeared. vignette("spectral") has the measured rates.
What it is not
The likelihood treats ordinates as independent, which is exact only
in the limit. It is biased for a short series or a spectrum near a
unit root; vignette("spectral") measures how much and says where it
is safe. Tapering reduces that bias and correlates the rows.
Examples
set.seed(1)
y <- as.numeric(arima.sim(list(ar = 0.6), 512))
pg <- frm_periodogram(y)
# a nonparametric log spectrum: the smoothing parameter is estimated
# by the same Laplace marginal likelihood as any other smooth
fit <- frm(bf(pgram ~ s(freq, k = 8)), family = whittle(), data = pg)
fixef(fit)
#> Estimate Est.Error Q2.5 Q97.5
#> Intercept 0.02596872 0.06279848 -0.09711403 0.1490515
#> sfreq_1 -4.72232641 1.48852879 -7.63978923 -1.8048636
# an averaged periodogram needs the shape it was averaged with
pg4 <- frm_periodogram(y, segments = 4)
frm(bf(pgram ~ log(freq)), family = whittle(tapers = 4), data = pg4)
#> Family: whittle
#> Links: mu = log; shape = log
#>
#> Formula: pgram ~ log(freq)
#> Data: pg4 (Number of observations: 63)
#> Method: ML logLik: -43.121 AIC: 90.242 BIC: 94.5283
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept -1.57 0.15 -1.86 -1.28 -10.62 <2e-16
#> logfreq -0.98 0.08 -1.14 -0.82 -12.15 <2e-16
#>
#> Fixed dpar: shape = 4