## Anomalous diffusion: Individual particle tracking vs. Ensemble averages

A recent article by Bos et al. (2017)^{1} in which the authors compared network measures derived from cross-sectional (between-individual) and within-individual data sparked some discussion on facebook.

The authors conclude that the outcomes of network analyses based on time-ordered ensemble averages are not the same as the outcomes based on the time evolution of individual trajectories.

I argued this result is to be expected, because I think the data-generating process underlying time and trial series of human physiology and performance should be considered a *non-ergodic process* and the particular analyses used to construct the graphical model (the Gaussian Graphical Model) assume ergodic processes. Simply put, the ergodic assumption is that the space-averaged behaviour of some variable observed in an ensemle will be the same as the time-averaged behaviour observed in an individual (… in the limit of infite time and space). If it is violated, this is a problem for such models, because it implies the statistical properties of the process change over time and one cannot depend (fullly) on universal laws of probability to predict or infer properties and behaviour of a system.

This post originated as a short tutorial on analysing stationarity and other properties of time series, the current text is an expansion of the page kindly tweeted by Eiko Fried:

Added to the tutorials section of psych-networks. If you have any other tutorials that would fit please let me knowhttps://t.co/zgJutWpOs8

— Eiko Fried (@EikoFried) May 19, 2017

### Particles are individuals too!

I often use the Complex Sysems argument to suggest measurements should be considered non-ergodic (many interactions between processes across multiple scales that are mostly multiplicative instead of additive), however, it turns out that studies of individual particle tracking (both biological and physical ‘particles’) also have to deal with *ergodicity-breaking*, *ageing* and *non-stationarity*. This behaviour is called *anomalous diffusion*, which turns out to be more common in nature than the classical *normal diffusion* processes studied by Brown, Einstein and many other prominent scientists.

The following quotes are from the article **Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking** which provides an excellent review and perspective.^{2} The authors list the conditions for normal diffusion, they should sound familiar because these assunmptions are the same as for most statistical analyses used in the social and life sciences:

The conditions assumed by Einstein in his derivation of the diffusion equation are (i) the independence of individual particles, (ii) the existence of a sufficiently small time scale beyond which individual displacements are statistically independent, and (iii) the property that the particle displacements during this time scale correspond to a typical mean free path distributed symmetrically in positive or negative directions.

So, we need:

*Independence*of indivdual measurement objects in a sample*Memorylessness*regarding any interactions between individuals and their environment. The after-effects of interactions should be short-lived and not affect behaviour in the long run (no long-range correlations).*Randomly distributed*deviations from central tendency (Gaussian error distribution).

What happens when we violate these assumptions?

In anomalous diffusion processes, at least one of these fundamental assumptions is violated, and the strong convergence to the Gaussian according to the central limit theorem broken. In particular, by departing from one or more of the assumptions (i)–(iii), we find that there exist many different generalisations of the Einstein–Smoluchowski diffusion picture.

Ok, but what does that mean?

The fact that we consider this large range of anomalous diffusion processes is the

non-universal nature of anomalous diffusionitself. Once we leave the realm of Brownian motion, we lose the confines of the central limit theorem forcing the processes to converge to the Gaussian behaviour predicted by Einstein.[…]

Quite commonly such analyses of time series from experiment or simulations are performed in terms of time averaged observables, in particular, the time averaged MSD [

Mean squared displecement = Mean squared ‘error’]. We point out that the physical interpretation of the obtained behaviour of such time averages in terms of thetypically available ensemble approaches may be treacherous: many of the anomalous diffusion processes discussed herein lead to adisparity between the ensemble and the time averaged observable, for instance, between the ensemble and time averaged MSDs […] even in the limit of long measurement times. Moreover, it turns out that individual results for time averages […] appear to beirreproducible, despite long measurement times.

To summarise, when non-ergodicity, non-stationarity and ageing play a role in the timeseries measurements of variables of interest, one cannot rely on ensemble statistics to yield valid inferences, in fact they should be considered **treacherous**, and individual time series averages are **irreproducible**.

What follows are some tests that have been developped (in different scientific disciplines) to check assumptions of non-stationarity (and more).

## A timeseries of daily mood: *down*

Unzip and load these interesting ESM data provided by:

Kossakowski, J. J., Groot, P. C., Haslbeck, J. M., Borsboom, D., & Wichers, M. (2016, December 12). Data from “Critical Slowing Down as a Personalized Early Warning Signal for Depression.” Retrieved from osf.io/j4fg8

There is a timecode in the dataset indicating when the participant was promted to nswer some questions. Here, we use it to create a time series object and plot it.

```
library(lubridate)
beep <- dmy_hms(paste0(df.ts$date," ",df.ts$beeptime,":00"),tz = "Europe/Amsterdam")
y <- xts(df.ts$mood_down, order.by = beep)
plot(y, main = "Time series of variable 'mood_down'")
```

## Test for **level** and **trend** stationarity

```
# Load some Time Series libraries
library(nonlinearTseries)
library(randtests)
library(tseries)
library(TSA)
```

```
# BARTELS RANK TEST [rank version of von NEUMANN's Ratio Test for Randomness]
# H0: randomness
# H1: nonrandomness
randtests::bartels.rank.test(df.ts$mood_down, alternative = "two.sided")
```

```
>
> Bartels Ratio Test
>
> data: df.ts$mood_down
> statistic = -16.606, n = 1474, p-value < 2.2e-16
> alternative hypothesis: nonrandomness
```

```
# H0: randomness
# H1: trend
randtests::bartels.rank.test(df.ts$mood_down, alternative = "left.sided")
```

```
>
> Bartels Ratio Test
>
> data: df.ts$mood_down
> statistic = -16.606, n = 1474, p-value < 2.2e-16
> alternative hypothesis: trend
```

```
# H0: randomness
# H1: systematic oscillation
randtests::bartels.rank.test(df.ts$mood_down, alternative = "right.sided")
```

```
>
> Bartels Ratio Test
>
> data: df.ts$mood_down
> statistic = -16.606, n = 1474, p-value = 1
> alternative hypothesis: systematic oscillation
```

```
# COX-STUART SIGN TEST
# H0: randomness
# H1: upward or downward trend
randtests::cox.stuart.test(na.exclude(df.ts$mood_down), alternative = "two.sided")
```

```
>
> Cox Stuart test
>
> data: na.exclude(df.ts$mood_down)
> statistic = 246, n = 362, p-value = 7.196e-12
> alternative hypothesis: non randomness
```

```
# H0: randomness
# H1: upward trend
randtests::cox.stuart.test(na.exclude(df.ts$mood_down), alternative = "right.sided")
```

```
>
> Cox Stuart test
>
> data: na.exclude(df.ts$mood_down)
> statistic = 246, n = 362, p-value = 3.598e-12
> alternative hypothesis: increasing trend
```

```
# H0: randomness
# H1: downward trend
randtests::cox.stuart.test(na.exclude(df.ts$mood_down), alternative = "left.sided")
```

```
>
> Cox Stuart test
>
> data: na.exclude(df.ts$mood_down)
> statistic = 246, n = 362, p-value = 1
> alternative hypothesis: decreasing trend
```

```
# KWIATKOWSKI-PHILLIPS-SCHMIDT-SHIN UNIT ROOT TEST ['reversed' DICKY-FULLER TEST]
# H0: level stationarity
# H1: unit root [not level stationary]
tseries::kpss.test(na.exclude(df.ts$mood_down), lshort = TRUE, null = "Level")
```

```
>
> KPSS Test for Level Stationarity
>
> data: na.exclude(df.ts$mood_down)
> KPSS Level = 1.7668, Truncation lag parameter = 7, p-value = 0.01
```

```
# H0: trend stationarity
# H1: not trend stationary
tseries::kpss.test(na.exclude(df.ts$mood_down), lshort = TRUE, null = "Trend")
```

```
>
> KPSS Test for Trend Stationarity
>
> data: na.exclude(df.ts$mood_down)
> KPSS Trend = 0.090975, Truncation lag parameter = 7, p-value = 0.1
```

## Test for **AR**, **ARCH** or an optimal **ARIMA** process

First inspect the partial autocorrelation function to get an idea of the AR order.

```
par(mfrow=c(1,2))
acf(df.ts$mood_down, na.action = na.pass)
pacf(df.ts$mood_down, na.action = na.pass)
```

One could conclude there’s at least \(AR(3)\) involved, however there are also possible long-range correlations in these data. I am not an initiate of the **how-to-interpret- (P)ACF-to-ARfiMA-orders-cult**, so let’s do some tests!

```
# KEENAN 1-DEGREE TEST OF NONLINEARITY
# H0: time series follows some AR process
# H1: time series cannot be considered some AR process
TSA::Keenan.test(na.exclude(df.ts$mood_down))
```

```
> $test.stat
> [1] 5.359887
>
> $p.value
> [1] 0.02074511
>
> $order
> [1] 16
```

Not an **AR** process, how about **ARCH** or a best fitting **ARiMA**?

```
# MCLEOD-LI TEST FOR CONDITIONAL HETEROSCEDASTICITY
# H0: time series follows some AR process
# H1: time series cannot be considered some ARCH process
TSA::McLeod.Li.test(y = df.ts$mood_down, plot = TRUE, omit.initial = TRUE)
```

```
# MCLEOD-LI TEST FOR CONDITIONAL HETEROSCEDASTICITY
# H0: time series follows some AR process
# H1: time series cannot be considered some ARiMA process
TSA::McLeod.Li.test(object = forecast::auto.arima(df.ts$mood_down), plot = TRUE, omit.initial = TRUE)
```

This test finds no convincing evidence of the timeseries being either **ARCH** or **ARiMA**.

## Structural decomposition

It is also possible to decompose the time series in different sources of variation. Let’s calculate the average frequency between beeps (in hours), and use that as a guess for any periodicity that may occur in the timeseries.

`(beepMean <- mean(time_length(diff(beep), unit = "hours")))`

`> [1] 3.886332`

```
ts0 <- ts(df.ts$mood_down, frequency = round(beepMean))
fit<-StructTS(ts0, type="BSM")
par(mfrow = c(4, 1)) # to give appropriate aspect ratio for next plot.
plot(cbind(fitted(fit), resids=resid(fit)), main= "Basic Structural Model for 'mood_down'")
```

There is a lot going on here, the residuals and level fluctuations show we probably should not assume a linear ergodic stochastic change process.

# Conclusions

This timeseries is:

- Not level stationary
- Not random
- Not oscillatory
- Not generated by some
**AR**process - Not
**ARCH**or generated by a best fitting**ARiMA**process - Contains long range correlations in
`PACF`

, e.g. lags 6, 9, 13, 25 and beyond (remember it is 25 * 4 hours!) - Upward Trend (which is stationary)

This is likely a violation of the strong ergodic condition / strong memorylessness property:

So, **do not** use any analytic techniques that depend on the ESM data and underlying process to be stationary and ergodic if these tests indicate the assumptions are violated.

**NEXT: Dynamics in state space**

What can we do to study non-ergodic processes?

Again from the anomalous diffussion article:

We note that non-ergodicity […] is not restricted to the spatial diffusion of particles, but similar principles hold for certain processes revealing non-exponential dynamics, such as the blinking behaviour of individual quantum dots or laser cooling. To physically interpret such measurements

we need to understand the time averages of individual time series. As we will see, this requires information beyond the conventional ensemble averages for a variety of anomalous diffusion processes.

This resonates with some great advice from Peter Molenaar in Todd Roses’s book The End of Average:

“first analyse, then aggregate”

When I find the time I’ll post some examples of how to represent and quantify state space dynamics of individual time series.

For now, the plots below should make clear the occurrence of state stransitions is not distributed according to a uniform or Gaussian pdf.

```
library(ggTimeSeries)
df.ts2 = data.frame(
Signal = round(head(df.ts$mood_down, -1),0),
NextSignal = round(tail(df.ts$mood_down, -1),0),
Weight = 1,
Label = "mood_down",
Size = 0,
Time = factor(round(seq(1,150, length.out = max(index(head(beep,-1))))))
)
N = length(df.ts2$Signal)
set.seed(123)
dfcat <- round(runif(N,min=-3,max=3))
dfcat2 <- data.frame(
Signal = round(head(dfcat, -1),0),
NextSignal = round(tail(dfcat, -1),0),
Weight = 1,
Label = "random_uniform",
Size=0,
Time = factor(round(seq(1,150, length.out = max(index(head(beep,-2))))))
)
dftot <- rbind(dfcat2,df.ts2)
ggplot(dftot, aes(xbucket = Signal,
ybucket = NextSignal,
fill = factor(NextSignal),
weight = Weight) )+
stat_marimekko(color = 'black', xlabelyposition = -0.1) +
scale_y_continuous("Relative frequency") +
scale_x_continuous("Current value") +
scale_fill_discrete("Next value") +
facet_wrap(~Label) +
ggtitle(label = "State tranisitions: Random numbers vs. 'I feel down'",subtitle = paste("box dimensions represent relative frequencies of",N,"values")) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
strip.background = element_rect(fill="grey90"),
strip.text = element_text(face = "bold"),
plot.background = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid = element_blank())
```

If we consider the values \(-3\) to \(3\) to span the possibe states the ‘system’ can be in, then a crude way to study the dynamics of the states is to look at a lag-1 *return plot*.

```
library(gganimate)
tmp <- as.data.frame(table(df.ts2$Signal,df.ts2$NextSignal))
# Set the size to frequency of occurrence of each combination of values
for(c in seq_along(tmp$Freq)){
ID <- (df.ts2$Signal%in%tmp$Var1[c])&(df.ts2$NextSignal%in%tmp$Var2[c])
if(sum(ID)>0){
df.ts2$Size[ID] <- tmp$Freq[c]
}
}
ga <- ggplot(df.ts2,aes(x=Signal,y=NextSignal, frame = Time)) +
geom_path(colour = "grey80") +
geom_point(aes(size=Size)) +
scale_y_continuous("Next value") +
scale_x_continuous("Current value") +
scale_size_continuous("Frequency of\n pair (Curren,Next)") +
ggtitle(label = "Window", subtitle = "Return plot") +
theme_bw() + coord_equal()
gganimate(ga, interval=.2, "returnplot.gif")
```

`#htmltools::HTML('{{< figure src = "returnplot.gif" >}}')`

It is appears to be the case the system is *attracted* to specific states. Perhaps the dynamics can be described as transitions between more or less stable states, or stable temporal patterns of recurring states.

to be continued …

天

Bos, F. M., Snippe, E., de Vos, S., Hartmann, J. A., Simons, C. J., van der Krieke, L., & Wichers, M. (2017). Can We Jump from Cross-Sectional to Dynamic Interpretations of Networks? Implications for the Network Perspective in Psychiatry.

*Psychotherapy and Psychosomatics, 86(3)*, 175-177↩Metzler, R., Jeon, J. H., Cherstvy, A. G., & Barkai, E. (2014). Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking.

*Physical Chemistry Chemical Physics, 16(44)*, 24128-24164.↩