Solved – How to perform a spectral density analysis in R using the multitaper package

rtime series

To Burr and other experts on this topic: I tested, well rather ad-hocly the R-package multitaper, and I succeeded getting graphical outputs. I got graphs for (a) spectrum, (b) spectrum with confidence interval and a graph for the F-test.
I have a couple of questions related to the package and the output which I hope will be commented.

  1. What is the F-test testing – what is the exact null-hypothesis (I suppose the null is the “white-noise” hypothesis – a flat spectral density)?
  2. Why is not the critical value for the F-test integrated in the F-test plot (as a horizontal line)? The package does not produce any numerical output from the test so we have no idea how “close” to reject the null-hypothesis.
  3. Is the significance level in the F-test optional, for example 5 or 10%?
  4. When you evaluate the spectrum, how do you transform frequencies to years when you have yearly data?
  5. It is also possible to plot the spectrum with confidence interval. What is the big deal of showing the confidence interval without relating it to the null-hypothesis?

  6. You recommended plotting the spectrum from spec.mtm on a log scale. Is it possible to integrate the log-command in the code: 

     resSpec <- spec.mtm(willamette, k=10, nw=5.0, nFFT = "default",
            centreWithSlepians = TRUE, Ftest = TRUE,
            jackknife = FALSE, maxAdaptiveIterations =100, 
            plot = TRUE, na.action = na.fail)

and how do we do it?

Best Answer

Just noticed this question. Here's responses to some of your questions.

The F-test test is (and I quote)

an F variance-ratio test with 2 and 2K-2 degrees of freedom for the significance of the estimated line component.

In practice, the test compares the value of the background spectra with the power in a line component, resulting in the above test. High values, high(er) significance.

As far as why the critical value is not included in the plot, I assume it is largely because which critical value you choose varies from problem to problem. The package does provide numerical output from the test (it is contained in the object returned from the function call, and you can pull it out and look at it as you wish). It even returns the two (numerator, denominator) degrees of freedom so you can compute the cutoff for any significance you want, as the distribution is essentially F, and has a known CDF.

To obtain the pieces you're looking for, pull the mtm object out of resSpec and take a look at it. For example:

spect<-resSpec$mtm
attributes(spect)
dof1<-2
dof2<-2*(spect$k)-2
Fval<-spect$Ftest

I'm not sure what you're asking when you say

Is the significance level in the F-test optional, for example 5 or 10%?

If you have data that is not sampled at $\delta t = 1s$, you can change the frequency delta by examining the spec class attached to resSpec. It is a time-series object, and can be modified in the usual way. Alternatively, if you know what you want, just create your own frequency vector and then plot(fr,resSpec$spec,etc).

The confidence interval that you can display on the spectrum is a jackknifed confidence interval. I can provide references if you'd like to read up on this a little more, but for most purposes, it doesn't provide directly applicable information for standard analysis.

As far as I know, there is no direct way to integrate the log spectrum command for the initial plot. However, I normally just turn auto-plot off and grab the resSpec$spec object and then manipulate it directly, fixing the $\delta f$ and setting the plot to plot(fr,log(resSpec$spec),etc).

Hope these answers helped.

Related Question