Fluctuation and Disperion analyses

Three plus two time series

We are going to analyse the same 3 time series as during the last session, now using fluctuation and dispersion analyses.

Importing data

Two ways:

A. By downloading:

  1. Follow the link, e.g. for series.xlsx.
  2. On the Github page, find a button marked Download (or Raw for text files).
  3. Download the file
  4. Load it into R using the code below

B. By directly importing the file in R from Github:

  1. Copy the url associated with the Download button on Github (right-click).
  2. The copied path should contain the word ‘raw’ somewhere in the url.
  3. Call rio::import(url)

Remember this is what the time series and retun plots looked like.

Now we’ll add two additional time series using the function noise_powerlaw(), run the code below:

4.0.1 Standardised Dispersion Analysis

Questions

  • Common data preparation before running fluctuation analyses like SDA:
    • Normalize the time series (using sd based on N, not N-1, e.g. by using ts_standardise(), or use the arguments of the function fd_sda())
    • Check whether time series length is a power of 2. If you use nextn(length(TS),factors=2), the number you get is the length your time series should be to create equal bin sizes. You can pad the series with 0s before analysis, look at package invctr at the infix function %+]%.
  • Perform sda on the 5 time series, use minData = 6 to exclude the larger bin sizes from estimating the scaling relation.
  • Compare what you find for fd_sda() to the other techniques (fd_RR,SampEn, acf, pacf)

Answers

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.64 | FD = 1.64 
## 
##  Fit range (n = 8)
## Slope = -0.56 | FD = 1.56
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdSDA2 <- fd_sda(series$TS_2, minData = 6, doPlot = TRUE, tsName = "TS2")
## 
## 
## fd_sda:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.58 | FD = 1.58 
## 
##  Fit range (n = 8)
## Slope = -0.39 | FD = 1.39
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdSDA3 <- fd_sda(series$TS_3, minData = 6, doPlot = TRUE, tsName = "TS3")
## 
## 
## fd_sda:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.64 | FD = 1.64 
## 
##  Fit range (n = 8)
## Slope = -0.59 | FD = 1.59
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdSDA4 <- fd_sda(series$TS_4, minData = 6, doPlot = TRUE, tsName = "TS4")
## 
## 
## fd_sda:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.75 | FD = 1.75 
## 
##  Fit range (n = 8)
## Slope = -0.18 | FD = 1.18
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdSDA5 <- fd_sda(series$TS_5, minData = 6, doPlot = TRUE, tsName = "TS5")
## 
## 
## fd_sda:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.05 | FD = 1.05 
## 
##  Fit range (n = 8)
## Slope = -0.03 | FD = 1.03
## 
## ~~~o~~o~~casnet~~o~~o~~~

# If you ask to return the powerlaw (returnPLAW = TRUE) the function plotFD_loglog() can make the plot as well
fdS5 <- fd_sda(series$TS_5, returnPLAW = TRUE, minData = 6, doPlot = TRUE, tsName = "TS5")
## 
## 
## fd_sda:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.05 | FD = 1.05 
## 
##  Fit range (n = 8)
## Slope = -0.03 | FD = 1.03
## 
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS5, title = "SDA", subtitle = "TS5")

4.0.2 Spectral Slope

Questions

  • Common data preparation before running spectral analyses:
    • Normalize the time series (using sd based on N, not N-1, e.g. by using ts_standardise(), or use the arguments of the function fd_psd())
    • Check whether time series length is a power of 2. If you use nextn(length(TS),factors=2), the number you get is the length your time series should be to create equal bin sizes. You can pad the series with 0s before analysis, look at package invctr at the infix function %+]%.
  • Perform fd_psd() on the 3 time series.
  • Compare to what you find for fd_psd to the other measures (fd_RR,SampEn, acf, pacf)

Answers

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = 0.18 | FD = 1.57 
## 
##  Hurvich-Deo (n = 36)
## Slope = 0.29 | FD = 1.61
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdPSD2 <- fd_psd(series$TS_2, doPlot = TRUE, tsName = "TS2")
## 
## 
## fd_psd:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = -0.27 | FD = 1.4 
## 
##  Hurvich-Deo (n = 79)
## Slope = -0.98 | FD = 1.2
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdPSD3 <- fd_psd(series$TS_3, doPlot = TRUE, tsName = "TS3")
## 
## 
## fd_psd:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = 0.12 | FD = 1.54 
## 
##  Hurvich-Deo (n = 32)
## Slope = 0.02 | FD = 1.51
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdPSD4 <- fd_psd(series$TS_4, doPlot = TRUE, tsName = "TS4")
## 
## 
## fd_psd:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = -0.84 | FD = 1.23 
## 
##  Hurvich-Deo (n = 28)
## Slope = 0.19 | FD = 1.57
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdPSD5 <- fd_psd(series$TS_5, doPlot = TRUE, tsName = "TS5")
## 
## 
## fd_psd:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = -1.9 | FD = 1.1 
## 
##  Hurvich-Deo (n = 51)
## Slope = -1.76 | FD = 1.11
## 
## ~~~o~~o~~casnet~~o~~o~~~

4.0.3 Detrended Fluctuation Analysis

Questions

  • Data preparation before running DFA analyses are not really necessary, except perhaps:
    • Check whether time series length is a power of 2. If you use nextn(length(TS),factors=2), the number you get is the length your time series should be to create equal bin sizes. You can pad the series with 0s before analysis, look at package invctr at the infix function %+]%.
  • Perform fd_dfa() on the 5 time series.
  • Compare to what you find for fd_dfa to the other measures (fd_RR,SampEn, acf, pacf)

Answers

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Detrended FLuctuation Analysis 
## 
##  Full range (n = 8)
## Slope = 0.44 | FD = 1.55 
## 
##  Exclude large bin sizes (n = 7)
## Slope = 0.45 | FD = 1.54
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdDFA2 <- fd_dfa(series$TS_2, doPlot = TRUE, tsName = "TS2")
## 
## 
## fd_dfa:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Detrended FLuctuation Analysis 
## 
##  Full range (n = 8)
## Slope = 1.1 | FD = 1.16 
## 
##  Exclude large bin sizes (n = 7)
## Slope = 1.23 | FD = 1.12
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdDFA3 <- fd_dfa(series$TS_3, doPlot = TRUE, tsName = "TS3")
## 
## 
## fd_dfa:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Detrended FLuctuation Analysis 
## 
##  Full range (n = 8)
## Slope = 0.47 | FD = 1.52 
## 
##  Exclude large bin sizes (n = 7)
## Slope = 0.47 | FD = 1.53
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdDFA4 <- fd_dfa(series$TS_4, doPlot = TRUE, tsName = "TS4")
## 
## 
## fd_dfa:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Detrended FLuctuation Analysis 
## 
##  Full range (n = 8)
## Slope = 0.95 | FD = 1.22 
## 
##  Exclude large bin sizes (n = 7)
## Slope = 0.98 | FD = 1.21
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdDFA5 <- fd_dfa(series$TS_5, doPlot = TRUE, tsName = "TS5")
## 
## 
## fd_dfa:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Detrended FLuctuation Analysis 
## 
##  Full range (n = 8)
## Slope = 1.35 | FD = 1.1 
## 
##  Exclude large bin sizes (n = 7)
## Slope = 1.4 | FD = 1.09
## 
## ~~~o~~o~~casnet~~o~~o~~~

4.0.4 Compare results

Questions

  • Create a table that compares the Fractal Dimension estimates based on the limited fit range for all time series (columns) with for each analysis (rows)
  • Create the same for the full fitrange.
  • What can you conclude about these series based on these results?

4.1 Surrogate testing: Beyond the straw-man null-hypothesis

Look at the explanation of surrogate analysis on the TiSEAN website (Here is direct link to the relevant sections)

4.1.1 Constrained realisations of complex signals

We’ll look mostly at three different kinds of surrogates:

  • Randomly shuffled: \(H_0:\) The time series data are independent random numbers drawn from some probability distribution.
  • Phase Randomised: \(H_0:\) The time series data were generated by a stationary linear stochastic process
  • Amplitude Adjusted ,Phase Randomised: \(H_0:\) The time series data were generated by a rescaled linear stochastic process.

Questions

  • Figure out how many surrogates you minimally need for a 2-sided test.
  • Package nonlinearTseries has a function FFTsurrogate() and package fractal has surrogate(). Look them up and try to create a surrogate test for time series TS1, TS3, and TS4 of the previous assignments.
    • If you use fractal::surrogate(), choose either aaft (amplitude adjusted Fourier transform) or phase (phase randomisation)
    • nonlinearTseries::FFTsurrogates will calculate phase randomised surrogates.
  • In package casnet there’s a function that will calculate and display a point-probability based on a distribution of surrogate measures and one observed measure: plotSUR_hist

Answers

Create a number of surrogates, calulate the slopes and plot the results.

  • First let’s get the values we found earlier (we’ll use the FD from sda, then we do the same for psd,but you can take any of the measures as well)
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.64 | FD = 1.64 
## 
##  Fit range (n = 9)
## Slope = -0.61 | FD = 1.61
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.64 | FD = 1.64 
## 
##  Fit range (n = 9)
## Slope = -0.58 | FD = 1.58
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.05 | FD = 1.05 
## 
##  Fit range (n = 9)
## Slope = -0.04 | FD = 1.04
## 
## ~~~o~~o~~casnet~~o~~o~~~
  • To create the surrogates using fractal::surrogate()
    • You’ll need to repeatedly call the function and store the result
    • To get into a dataframe we have to use unclass() otherwise R doesn’t recognise the object
    • Function t() transposes the data from 39 rows to 39 columns
  • Once you have the surrogates in a datafram, repeatedly calculate the measure you want to compare, in this case fd_sda().

  • Collect the results in a dataframe, so we can compare.
  • Call function plotSUR_hist() to see the results of the test.

  • Let’s print the results into a table
Source FDsda FDsdaAAFT.median
TS1 1.61 1.59
TS3 1.58 1.55
TS5 1.04 1.07

FD: SPECTRAL SLOPE

## 
## 
## fd_psd:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = 0.18 | FD = 1.57 
## 
##  Hurvich-Deo (n = 36)
## Slope = 0.29 | FD = 1.61
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
## 
## fd_psd:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = 0.12 | FD = 1.54 
## 
##  Hurvich-Deo (n = 32)
## Slope = 0.02 | FD = 1.51
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
## 
## fd_psd:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = -1.91 | FD = 1.1 
## 
##  Hurvich-Deo (n = 73)
## Slope = -1.92 | FD = 1.1
## 
## ~~~o~~o~~casnet~~o~~o~~~
  • To create the surrogates using fractal::surrogate()
    • You’ll need to repeatedly call the function and store the result
    • To get into a dataframe we have to use unclass() otherwise R doesn’t recognise the object
    • Function t() transposes the data from 39 rows to 39 columns
  • Once you have the surrogates in a dataframe, repeatedly calculate the measure you want to compare, in this case fd_sda().
  • Call function plotSUR_hist() to get the results of the test.

  • Let’s print the results into a table
Source FDsda FDsdaAAFT.median FDpsd FDpsdAAFT.median
TS1 1.61 1.59 1.61 1.59
TS3 1.58 1.55 1.51 1.55
TS5 1.04 1.07 1.10 1.13