Fluctuation and Disperion analyses

Importing data

Two ways:

A. By downloading:

  1. Follow the link, e.g. for noises.xlsx.
  2. On the Github page, find a button marked Download (or Raw for textfiles).
  3. Download the file
  4. Load it into R using the code below
library(rio)
noises <- rio::import("noises.xlsx")

B. By importing from Github:

  1. Copy the url associated with the Download button on Github (right-clik).
  2. The copied path should contain the word ‘raw’ somewhere in the url.
  3. Call rio::import(url)
library(rio)
noises <- rio::import("https://github.com/complexity-methods/CSA-assignments/raw/master/assignment_data/BasicTSA_arma/noises.xlsx")

2.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 function arguments)
    • Check whether time series length is a power of 2. If you use N <- log2(length(y)), the number you get is 2^N. You need an integer power for the length in order to create equal bin sizes. You can pad the series with 0s, or make it shorter before analysis.
  • Perform sda on the 3 time series, or any of the other series you already analysed.
  • Compare to what you find for fd_sda to the other techniques (fd_RR,SampEn, acf, pacf)

Answers

library(casnet)

fdS1 <- fd_sda(noises$TS_1, returnPLAW = TRUE)
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.4 | FD = 1.4 
## 
##  Fit range (n = 9)
## Slope = -0.15 | FD = 1.15
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdS2 <- fd_sda(noises$TS_2, returnPLAW = TRUE)
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.54 | FD = 1.54 
## 
##  Fit range (n = 9)
## Slope = -0.5 | FD = 1.5
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdS3 <- fd_sda(noises$TS_3, returnPLAW = TRUE)
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -1.39 | FD = 2.39 
## 
##  Fit range (n = 9)
## Slope = -1.46 | FD = 2.46
## 
## ~~~o~~o~~casnet~~o~~o~~~

# If you asked to return the powerlaw (returnPLAW = TRUE) the function plotFD_loglog() can make a plot
plotFD_loglog(fdS1,title = "SDA", subtitle = "Series 1")

plotFD_loglog(fdS2,title = "SDA", subtitle = "Series 2")

plotFD_loglog(fdS3,title = "SDA", subtitle = "Series 3")

2.0.2 Spectral Slope

Questions

  • Common data preparation before running spectral analyses:
    • Normalize the time series (using sd based on N, not N-1)
    • Check whether time series length is a power of 2. If you use N <- log2(length(y)), the number you get is 2^N. You need an integer power for the length in order to create equal bin sizes. You can pad the series with 0s, or make it shorter before analysis.
  • Perform fd_psd() on the 3 time series, or any of the other series you already analysed.
  • Compare to what you find for fd_psd to the other measures (fd_RR,SampEn, acf, pacf)

Answers


fdS1 <- fd_psd(noises$TS_1, doPlot = TRUE)
## 
## 
## fd_psd:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = -1.07 | FD = 1.19 
## 
##  Hurvich-Deo (n = 61)
## Slope = -1.24 | FD = 1.16
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdS2 <- fd_psd(noises$TS_2, returnPLAW = TRUE)
## 
## 
## fd_psd:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = 0 | FD = 1.5 
## 
##  Hurvich-Deo (n = 54)
## Slope = -0.06 | FD = 1.48
## 
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS2,subtitle = "Series 2")

fdS3 <- fd_psd(noises$TS_3, returnPLAW = TRUE)
## 
## 
## fd_psd:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = 2.02 | FD = 1.9 
## 
##  Hurvich-Deo (n = 21)
## Slope = -2.15 | FD = 1.09
## 
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS3,subtitle = "Series 3")

2.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 N <- log2(length(y)), the number you get is 2^N. You need an integer power for the length in order to create equal bin sizes. You can pad the series with 0s, or make it shorter before analysis.
  • Perform fd_dfa() on the 3 time series, or any of the other series you already analysed.
  • Compare to what you find for fd_dfa to the other measures (fd_RR,SampEn, acf, pacf)

Answers


fdS1 <- fd_dfa(noises$TS_1, doPlot = TRUE)
## 
## 
## fd_dfa:  Sample rate was set to 1.

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Detrended FLuctuation Analysis 
## 
##  Full range (n = 8)
## Slope = 1.02 | FD = 1.19 
## 
##  Exclude large bin sizes (n = 7)
## Slope = 1.02 | FD = 1.19
## 
## ~~~o~~o~~casnet~~o~~o~~~
fdS2 <- fd_dfa(noises$TS_2, returnPLAW = TRUE)
## 
## 
## fd_dfa:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Detrended FLuctuation Analysis 
## 
##  Full range (n = 8)
## Slope = 0.49 | FD = 1.51 
## 
##  Exclude large bin sizes (n = 7)
## Slope = 0.52 | FD = 1.48
## 
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS2,subtitle = "Series 2")

fdS3 <- fd_dfa(noises$TS_3, returnPLAW = TRUE)
## 
## 
## fd_dfa:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Detrended FLuctuation Analysis 
## 
##  Full range (n = 8)
## Slope = 0.02 | FD = 1.98 
## 
##  Exclude large bin sizes (n = 7)
## Slope = 0.02 | FD = 1.98
## 
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS3,subtitle = "Series 3")

2.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)

2.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 surrogate. Look them up and try to create a surrogate test for some of the time series 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)
library(fractal)

# RR assignment data `TS1,TS2,TS3`
# Now we'll trim the time noises length to a power of 2 to get the cleanest results, this is not strictly necessary 
L1 <- floor(log2(length(noises$TS_1)))
L2 <- floor(log2(length(noises$TS_2)))
L3 <- floor(log2(length(noises$TS_3)))

TS1_FDsda <- fd_sda(noises$TS_1[1:2^L1])$fitRange$FD
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.4 | FD = 1.4 
## 
##  Fit range (n = 9)
## Slope = -0.15 | FD = 1.15
## 
## ~~~o~~o~~casnet~~o~~o~~~
TS2_FDsda <- fd_sda(noises$TS_2[1:2^L2])$fitRange$FD
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -0.54 | FD = 1.54 
## 
##  Fit range (n = 9)
## Slope = -0.5 | FD = 1.5
## 
## ~~~o~~o~~casnet~~o~~o~~~
TS3_FDsda <- fd_sda(noises$TS_3[1:2^L3])$fitRange$FD
## 
## 
## fd_sda:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Standardised Dispersion Analysis 
## 
##  Full range (n = 10)
## Slope = -1.39 | FD = 2.39 
## 
##  Fit range (n = 9)
## Slope = -1.46 | FD = 2.46
## 
## ~~~o~~o~~casnet~~o~~o~~~
# y1<-growth_ac(r = 2.9,type="logistic",N = 2^L1)
# TSlogMapr29 <- fd_sda(y1)$fitRange$FD
# 
# y2<-growth_ac(r = 4,type="logistic", N = 2^L1)
# TSlogMapr4 <- fd_sda(y2)$fitRange$FD
# 
# # These were rendomised
# TS1R_FDsda <- fd_sda(TS1Random[1:2^L1])$fitRange$FD
# TS2R_FDsda <- fd_sda(TS1Random[1:2^L2])$fitRange$FD
# TS3R_FDsda <- fd_sda(TS1Random[1:2^L3])$fitRange$FD
  • 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() trnasposes 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().
# NOW CREATE SURROGATES
library(plyr)
library(tidyverse)
library(fractal)

# For a two-sided test at alpha = .05 we need N=39
Nsurrogates <- 39

TS1surrogates <- as.data.frame(t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=noises$TS_1[1:2^L1] ,method="aaft")))))
colnames(TS1surrogates) <- paste0("S",1:NCOL(TS1surrogates))

# Add the observed data
TS1surrogates$Obs <- noises$TS_1[1:2^L1]
plotTS_multi(TS1surrogates)

library(fractal)
# Now we calculate FD for each series
TS1surrogates_FD <- plyr::laply(1:Nsurrogates,function(s) fd_sda(y=TS1surrogates[,s], silent = TRUE)$fitRange$FD)

TS2surrogates <- t(plyr::ldply(1:Nsurrogates, function(s) unclass(surrogate(x=noises$TS_2[1:2^L2] ,method="aaft"))))
colnames(TS2surrogates) <- paste0("aaftSurr",1:NCOL(TS2surrogates))

TS2surrogates_FD <- plyr::laply(1:Nsurrogates,function(s) fd_sda(y=TS2surrogates[,s], silent = TRUE)$fitRange$FD)


TS3surrogates <- t(plyr::ldply(1:Nsurrogates, function(s) unclass(surrogate(x=noises$TS_3[1:2^L3] ,method="aaft"))))
colnames(TS3surrogates) <- paste0("aaftSurr",1:NCOL(TS3surrogates))

TS3surrogates_FD <- plyr::laply(1:Nsurrogates,function(s) fd_sda(y=TS3surrogates[,s], silent = TRUE)$fitRange$FD)

# 
# TSr29surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=y1, method="aaft"))))
# colnames(TSr29surrogates) <- paste0("aaftSurr",1:NCOL(TSr29surrogates))
# 
# TSr29surrogates_FD <- laply(1:Nsurrogates,function(s) fd_sda(y=TSr29surrogates[,s], silent = TRUE)$fitRange$FD)
# 
# 
# TSr4surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=y2, method="aaft"))))
# colnames(TSr4surrogates) <- paste0("aaftSurr",1:NCOL(TSr4surrogates))
# 
# TSr4surrogates_FD <- laply(1:Nsurrogates,function(s) fd_sda(y=TSr4surrogates[,s], silent = TRUE)$fitRange$FD)
  • Collect the results in a dataframe, so we can compare.
x <- data.frame(Source = c("TS1", "TS2", "TS3"), # "TSlogMapr29", "TSlogMapr4"), 
                FDsda = c(TS1_FDsda, TS2_FDsda, TS3_FDsda), # TSlogMapr29, TSlogMapr4), 
                #FDsdaRandom =  c(TS1R_FDsda, TS2R_FDsda, TS3R_FDsda,NA,NA), 
                FDsdaAAFT.median= c(median(TS1surrogates_FD),median(TS2surrogates_FD),median(TS3surrogates_FD))) #,median(TSr29surrogates_FD),median(TSr4surrogates_FD)))
  • Call function plotSUR_hist() to see the results of the test.
plotSUR_hist(surrogateValues = TS1surrogates_FD, observedValue = TS1_FDsda, sides = "two.sided", doPlot = TRUE, measureName = "sda FD TS1")

plotSUR_hist(surrogateValues = TS2surrogates_FD, observedValue = TS2_FDsda, sides = "two.sided", doPlot = TRUE, measureName = "sda FD TS2")

plotSUR_hist(surrogateValues = TS3surrogates_FD, observedValue = TS3_FDsda,sides = "two.sided", doPlot = TRUE, measureName = "sda FD TS3")

# 
# plotSUR_hist(surrogateValues = TSr29surrogates_FD, observedValue = TSlogMapr29,sides = "two.sided", doPlot = TRUE, measureName = "sda FD logMap r=2.9")
# 
# plotSUR_hist(surrogateValues = TSr4surrogates_FD, observedValue = TSlogMapr4, sides = "two.sided", doPlot = TRUE, measureName = "sda FD logMap r=4")

FD: SPECTRAL SLOPE

# RR assignment data `TS1,TS2,TS3`
L1 <- floor(log2(length(noises$TS_1)))
L2 <- floor(log2(length(noises$TS_2)))
L3 <- floor(log2(length(noises$TS_2)))

TS1_FDpsd <- fd_psd(noises$TS_1[1:2^L1])$fitRange$FD
## 
## 
## fd_psd:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = -1.07 | FD = 1.19 
## 
##  Hurvich-Deo (n = 61)
## Slope = -1.24 | FD = 1.16
## 
## ~~~o~~o~~casnet~~o~~o~~~
TS2_FDpsd <- fd_psd(noises$TS_2[1:2^L2])$fitRange$FD
## 
## 
## fd_psd:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = 0 | FD = 1.5 
## 
##  Hurvich-Deo (n = 54)
## Slope = -0.06 | FD = 1.48
## 
## ~~~o~~o~~casnet~~o~~o~~~
TS3_FDpsd <- fd_psd(noises$TS_3[1:2^L3])$fitRange$FD
## 
## 
## fd_psd:  Sample rate was set to 1.
## 
## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
##  Power Spectral Density Slope 
## 
##  All frequencies (n = 512)
## Slope = 2.02 | FD = 1.9 
## 
##  Hurvich-Deo (n = 21)
## Slope = -2.15 | FD = 1.09
## 
## ~~~o~~o~~casnet~~o~~o~~~
# # Logistic map - make it the length of L1
# y1<-growth_ac(r = 2.9,type="logistic",N = 2^L1)
# # Turn standardisation and detrending off!
# TSlogMapr29 <- fd_psd(y1, standardise = FALSE, detrend = FALSE)$fitRange$FD
# y2<-growth_ac(r = 4,type="logistic", N = 2^L1)
# TSlogMapr4 <- fd_psd(y2, standardise = FALSE, detrend = FALSE)$fitRange$FD
# 
# TS1R_FDpsd <- fd_psd(TS1Random[1:2^L1])$fitRange$FD
# TS2R_FDpsd <- fd_psd(TS1Random[1:2^L2])$fitRange$FD
# TS3R_FDpsd <- fd_psd(TS1Random[1:2^L3])$fitRange$FD
  • 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() trnasposes 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().
# NOW CREATE SURROGATES
library(dplyr)

# For a two-sided test at alpha = .05 we need N=39
Nsurrogates <- 39

TS1surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(fractal::surrogate(x=noises$TS_1[1:2^L1] ,method="aaft"))))
colnames(TS1surrogates) <- paste0("aaftSurr",1:NCOL(TS1surrogates))

TS1surrogates_FD <- laply(1:Nsurrogates,function(s){fd_psd(y=TS1surrogates[,s], silent = TRUE)$fitRange$FD})


TS2surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=noises$TS_2[1:2^L2] ,method="aaft"))))
colnames(TS2surrogates) <- paste0("aaftSurr",1:NCOL(TS2surrogates))

TS2surrogates_FD <- laply(1:Nsurrogates,function(s) fd_psd(y=TS2surrogates[,s], silent = TRUE)$fitRange$FD)


TS3surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=noises$TS_3[1:2^L3] ,method="aaft"))))
colnames(TS3surrogates) <- paste0("aaftSurr",1:NCOL(TS3surrogates))

TS3surrogates_FD <- laply(1:Nsurrogates,function(s) fd_psd(y=TS3surrogates[,s], silent = TRUE)$fitRange$FD)


# TSr29surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=y1, method="aaft"))))
# colnames(TSr29surrogates) <- paste0("aaftSurr",1:NCOL(TSr29surrogates))
# 
# TSr29surrogates_FD <- laply(1:Nsurrogates,function(s) fd_psd(y=TSr29surrogates[,s], silent = TRUE)$fitRange$FD)
# 
# 
# TSr4surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=y2, method="aaft"))))
# colnames(TSr4surrogates) <- paste0("aaftSurr",1:NCOL(TSr4surrogates))
# 
# TSr4surrogates_FD <- laply(1:Nsurrogates,function(s) fd_psd(y=TSr4surrogates[,s], silent = TRUE)$fitRange$FD)
  • Call function plotSUR_hist() to get the results of the test.
plotSUR_hist(surrogateValues = TS1surrogates_FD, observedValue = TS1_FDpsd, sides = "two.sided", doPlot = TRUE, measureName = "psd FD TS1")

plotSUR_hist(surrogateValues = TS2surrogates_FD, observedValue = TS2_FDpsd,sides = "two.sided", doPlot = TRUE, measureName = "psd FD TS2")

plotSUR_hist(surrogateValues = TS3surrogates_FD, observedValue = TS3_FDpsd,sides = "two.sided", doPlot = TRUE, measureName = "psd FD TS3")

# plotSUR_hist(surrogateValues = TSr29surrogates_FD, observedValue = TSlogMapr29,sides = "two.sided", doPlot = TRUE, measureName = "psd FD logMap r=2.9")
# 
# plotSUR_hist(surrogateValues = TSr4surrogates_FD, observedValue = TSlogMapr4, sides = "two.sided", doPlot = TRUE, measureName = "psd FD logMap r=4")
  • Let’s print the results into a table
x$FDpsd <- c(TS1_FDpsd,TS2_FDpsd,TS3_FDpsd) #TSlogMapr29,TSlogMapr4)
#x$FDpsdRandom <- c(TS1R_FDpsd, TS2R_FDpsd, TS3R_FDpsd,NA,NA)
x$FDpsdAAFT.median <- c(median(TS1surrogates_FD),median(TS2surrogates_FD),median(TS3surrogates_FD)) #,median(TSr29surrogates_FD),median(TSr4surrogates_FD))
knitr::kable(x, digits = 2, booktabs=TRUE,formt="html")
Source FDsda FDsdaAAFT.median FDpsd FDpsdAAFT.median
TS1 1.61 1.57 1.61 1.58
TS2 1.55 1.54 1.20 1.23
TS3 1.58 1.57 1.51 1.54