The fractalRegression package *cannot* currently be installed
through CRAN (the Comprehensive R Archive Network) https://cran.r-project.org/package=fractalRegression.
But, once it’s available there, it can be installed withe the
following:

`install.packages("fractalRegression")`

The development version is available on Aaron Likens’ Github (https://github.com/aaronlikens/fractalRegression) and can be installed using R devtools. This is a source package and requires compilation of C++ code. Windows users can install RTools software package to get necessary components: https://cran.r-project.org/bin/windows/Rtools/ Intel Mac users can install xcode along with the xcode commandline tools. Users with Mac silicon will need to a bit more fiddling to get the build chain working including the R recommended gfortran compiler. Some useful tips can be found here: https://mac.r-project.org/.

`::install_github("aaronlikens/fractalRegression") devtools`

The aim with this `fractalRegression`

package, and the
subsequent vignettes, is to show users how to use the various functions
in the package to perform univariate mono- and multi- fractal analysis
as well a bivariate fractal correlation and regression analyses.

In this first example, we apply DFA, originally developed by Peng et
al. (1994) to simulated white noise, pink noise, and anti-correlated
noise using the `dfa()`

function.

```
# Generate white noise
set.seed(23422345)
<- rnorm(5000)
noise <- logscale(scale_min=16, scale_max = 1025, scale_ratio = 2)
scales
# Run DFA on white noise
<- dfa(x = noise, order = 1, verbose = 1, scales=scales, scale_ratio = 2) dfa.noise.out
```

Since we simulated this white noise, we would expect that our
*α* is close to 0.5. Indeed, we observe the estimate from the
above analysis is 0.5207799. Next we use the `fgn_sim()`

to
simulate fractional Gaussian noise with known properties.

```
# Generate pink noise
set.seed(34534534)
<- fgn_sim(n = 5000, H = 0.9)
pink.noise
# Run DFA on pink noise
<- dfa(x = pink.noise, order = 1, verbose = 1, scales = scales, scale_ratio = 2) dfa.pink.out
```

Since we simulated this data with an H = 0.9, we would expect that
our *α* is close to this value. Again, we observe the estimate
from the above analysis is 0.833253. Next, we use `fgn_sim()`

to generate some anti-correlated noise.

```
# Generate anti-correlated noise
set.seed(24633453)
<- fgn_sim(n = 5000, H = 0.25)
anticorr.noise
# Run DFA on anti-correlated noise
<- dfa(x = anticorr.noise, order = 1, verbose = 1, scales = scales, scale_ratio = 2) dfa.anticorr.out
```

As with above, since we simulated with H = 0.25, we observed a close
estimate of the *α* exponent, which is 0.2301879.

Now let’s take a look at the log-log plots for the three types of
noise using the `dfa.plot()`

function. Given the estimates
above, we see that the slopes for white noise, pink noise, and
anti-correlated noise conform to our expectations. In the figures below,
the intercept and the slopes (i.e., *α* exponents) are shown in
addition to the *R*^{2}.

```
par(mfrow=c(3,1))
dfa.plot(dfa.noise.out)
dfa.plot(dfa.pink.out)
dfa.plot(dfa.anticorr.out)
```

Now that we have run a mono-fractal analysis using
`dfa()`

, we want to expand this to look for evidence of
multifractality using multifractal detrended fluctuation analysis
(MF-DFA) developed by Kantelhardt et al. (2002). That is, we aim to
determine whether there is evidence of multiple scaling relationships
and interactions across scales. We can do this easily using the
`mfdfa()`

function.

```
# Run MF-DFA on simulated pink and white noise
<- mfdfa(x = pink.noise, q = c(-5:5), order = 1, scales = scales)
mf.dfa.pink.out
<- mfdfa(x = noise, q = c(-5:5), order = 1, scale = scales) mf.dfa.white.out
```

Let’s first make sure that our *α* estimate is the same as
before when q = 2. We can check this easily with the code below, and it
looks good. For example, the pink noise when q = 2, Hq = 0.833253, which
is equal to our *α* = 0.833253 from above.

`$Hq[mf.dfa.pink.out$q==2] mf.dfa.pink.out`

`## [1] 0.833253`

`$Hq[mf.dfa.white.out$q==2] mf.dfa.white.out`

`## [1] 0.5207799`

Next we are going to work with data that we include in our package
(`fractaldata`

). This data was originally provided by Ihlen
(2012). It includes a white noise time series, monofractal time series,
and a multifractal time series.

Now let’s run MFDFA on these times series. In this case we replicate the choice of parameters in Ihlen (2012) with a q range of -5 to 5, and order = 1, a scale min of 16, and a scale max 1,024.

```
<- logscale(scale_min = 16, scale_max = 1025, scale_ratio = 1.1)
scales <- mfdfa(x = fractaldata$whitenoise, q = c(-5:5), order = 1, scales = scales)
white.mf.dfa.out
<- mfdfa(x = fractaldata$monofractal, q = c(-5:5), order = 1, scales = scales)
mono.mf.dfa.out
<- mfdfa(x = fractaldata$multifractal, q = c(-5:5), order = 1, scales = scales) multi.mf.dfa.out
```

A common way to understand if there is evidence of multifractality is
to examine a plot showing the Hq estimates for different values of q. If
all the plots have the same slope, that provides evidence of
monofractality. If there are distinct slopes, then there is evidence of
multifractality. It’s also important to check here that the slopes of
`log_2_scale`

and `log_fq`

are linear, thus
implying that they are scale invariant. If not, then it could be the
case that a higher order polynomial detrending is appropriate (see
Kantelhardt et al., 2001). We are going to use the
`mfdfa.plot()`

function to compare the monofractal and
multifractal signals.

```
# Let's first make a plot of the monofractal case
mfdfa.plot(mono.mf.dfa.out)
```

`## Loading required package: colorRamps`

Now let’s compare the above plot for the monofractal and multifractal results. In comparing, the top left part of both plots, we see that the slopes of the lines for the multifractal signal are divergent compared to the monofractal case.

`mfdfa.plot(multi.mf.dfa.out)`

It’s also common to examine the relationship between Hq and q as well as H and D(H). We can see the comparison in the bottom right of the two figures above, and the relative difference in the widths of the mutlifractal spectra.

A common metric for comparing the multifractal spectrum is to
calculate the width (W) as the
*h*_{max} − *h*_{min}.
Let’s do this to compare the monofractal and multifractal time series.
While from the figure above, it’s quite clear that the width of the
multifractal spectrum for the multifractal signal is larger, let’s
calculate it here anyways for the sake of completeness and because
spectrum width can be used as dependent variable in statistical
analyses.

```
<- max(mono.mf.dfa.out$h) - min(mono.mf.dfa.out$h)
mono.W <- max(multi.mf.dfa.out$h) - min(multi.mf.dfa.out$h) multi.W
```

We observe here that the W for the multifractal signal is 1.3566793 and for the monofractal signal, W is 0.0754753.

Moving on from the univariate analyses, if we have two non-stationary
signals and we want to examine the scale-wise fluctuations as well as
the scale-wise cross-correlation of these fluctuations, we can use DCCA
using the `dcca()`

function, which was originally developed
originally by Podobnik and Stanley (2008) and Zebende (2011).

First, however, we are going to simulate some data first. These
functions are taken from and available at Ladislav Kristoufek’s website
(http://staff.utia.cas.cz/kristoufek/Ladislav_Kristoufek/Codes_files/MC-ARFIMA.R)
and they correspond with the following paper (Kristoufek, 2013). We
implement a function that includes all of these options called
`mc_ARFIMA()`

.

The MC-ARFIMA models take the form of the two equations shown below:

\(x_t = \\alpha \\sum^{+ \\infty}\_{n=0}a_n(d_1)\\varepsilon\_{1,t-n}+\\beta \\sum^{+ \\infty}\_{n=0}a_n(d_2)\\varepsilon\_{2,t-n}\)

\(y_t = \\gamma \\sum^{+ \\infty}\_{n=0}a_n(d_3)\\varepsilon\_{3,t-n}+\\delta \\sum^{+ \\infty}\_{n=0}a_n(d_4)\\varepsilon\_{4,t-n}\)

In this case, we use the parameters from Kristoufec (2013) for Model 1 (p. 6487) in which case the resulting two time series of length 10,000 exhibit long range correlations (LRC) as well as long range cross-correlations (LRCC).

```
set.seed(987345757)
<- mc_ARFIMA(process='Mixed_ARFIMA_ARFIMA',alpha = 0.2, beta = 1, gamma = 1, delta = 0.2, n = 10000, d1 = 0.4, d2 = 0.3, d3 = 0.3, d4=0.4, rho=0.9)
sim1 plot(sim1[,1],type='l', ylab= "Signal Amplitude", xlab='Time', main = "MC-ARFIMA with LRC and LRCC")
lines(sim1[,2], col='blue')
```

```
<- logscale(scale_min = 10, scale_max = 1000, scale_ratio = 1.1)
scales <- dcca(sim1[,1], sim1[,2], order = 1, scales = scales)
dcca.out.arfima <- as.data.frame(dcca.out.arfima)
dcca.out.arfima <- sd(dcca.out.arfima$rho)/sqrt(length(dcca.out.arfima$rho))
error <- ggplot(data=dcca.out.arfima, aes(x=scales,y=rho)) + geom_point() +geom_line() + ggtitle("DCCA on MC-ARFIMA processes with LRC and LRCC")+ geom_pointrange(aes(ymin=rho-error,ymax=rho+error))
dcca.plot #geom_smooth(method=lm, formula = y ~ poly(x, 2), se = TRUE)
dcca.plot
```

In the above figure, we see that the correlation between the MC-ARFIMA processes are consistently high and continue to be high at increasing time scales. Standard errors are plotted around each point. Now let’s contrast this with MC-ARFIMA processes with LRC and short-range cross-correlation (SRCC).

```
set.seed(423422)
<- mc_ARFIMA(process = "Mixed_ARFIMA_AR", alpha = 1,beta = 1,gamma = 1,delta = 1,n =10000,d1=0.4,d2=0.4,theta1=0.8,theta2=0.8,rho=0.9)
sim2 plot(sim2[,1],type='l', ylab= "Signal Amplitude", xlab='Time', main = "MC-ARFIMA with LRC and SRCC")
lines(sim2[,2], col='blue')
```

```
<- logscale(scale_min = 10, scale_max = 1001, scale_ratio = 1.1)
scales <- dcca(sim2[,1], sim2[,2], order = 1, scales = scales)
dcca.out.arfima2 <- as.data.frame(dcca.out.arfima2)
dcca.out.arfima2 <- sd(dcca.out.arfima2$rho)/sqrt(length(dcca.out.arfima2$rho))
error <- ggplot(data=dcca.out.arfima2, aes(x=scales,y=rho)) + geom_point() +geom_line() + ggtitle("DCCA on MC-ARFIMA processes with LRC and SRCC") + geom_pointrange(aes(ymin=rho-error,ymax=rho+error))
dcca.plot2 #geom_smooth(method=lm, formula = y ~ poly(x, 2), se = TRUE)
dcca.plot2
```

In contrast to the previous DCCA analysis, the above figure shows a signal that begins with a high cross-correlation, but that begins to deviate and trend lower substantially at increasing scale sizes.

Let’s consider the above simulations, and consider the question: what
if we want to from a scale-wise correlation framework to a regression
framework? Or, put differently, can we use the scale-wise fluctuations
of one variable to predict the scale-wise fluctuations of the other? As
with a traditional regression approach, we will use one of our variables
as our predictor (*x*_{t}) and the other as our
outcome (*y*_{t}).

```
<- logscale(scale_min = 10, scale_max = 1001, scale_ratio = 1.1)
scales <- as.data.frame(mra(x = sim1[,1], y = sim1[,2],order = 1, scales = scales))
mra.out <- sd(mra.out$betas)/sqrt(length(mra.out$betas))
error <- ggplot(data=mra.out, aes(x=scales,y=betas)) + geom_point() +geom_line() +ggtitle("Multiscale Regression Analysis for MC-ARFIMA with LRC and LRCC") + geom_pointrange(aes(ymin=betas-error,ymax=betas+error))
mra.plot #geom_smooth(method=lm, formula = y ~ poly(x, 2), se = TRUE)
mra.plot
```

Generally, we observe that the *β* coefficients are relatively
stable at increasing time scales with a general, perhaps quadratically
increasing trend. Here it is also important to investigate the change in
*R*^{2} as well as the *t*-values. Below we see
that the *R*^{2} is quite high at most of the time scales
with *R*_{min}^{2}=
0.67 and *R*_{max}^{2}=
1.85and all of the *t*-values greater than the conventional
cut-off at 1.96. So between these two component ARFIMA processes, the
output of MRA shows that much of the scale specific variance in
*y*_{t} is explained and predicted by
*x*_{t}.

`::kable(mra.out, format='html', digits =2,align='c',caption = "Output from MRA") knitr`

scales | betas | r2 | t_observed |
---|---|---|---|

10 | 0.83 | 1.85 | NaN |

11 | 0.81 | 1.72 | NaN |

12 | 0.81 | 1.60 | NaN |

13 | 0.79 | 1.49 | NaN |

14 | 0.80 | 1.41 | NaN |

16 | 0.78 | 1.29 | NaN |

17 | 0.77 | 1.24 | NaN |

19 | 0.75 | 1.17 | NaN |

21 | 0.76 | 1.11 | NaN |

23 | 0.74 | 1.07 | NaN |

25 | 0.74 | 1.02 | NaN |

28 | 0.73 | 0.98 | 29.88 |

31 | 0.73 | 0.93 | 17.94 |

34 | 0.72 | 0.90 | 15.67 |

37 | 0.69 | 0.89 | 15.03 |

41 | 0.69 | 0.88 | 15.64 |

45 | 0.68 | 0.85 | 14.63 |

50 | 0.67 | 0.83 | 14.66 |

55 | 0.65 | 0.79 | 13.35 |

61 | 0.68 | 0.82 | 15.90 |

67 | 0.66 | 0.81 | 15.86 |

74 | 0.65 | 0.80 | 16.35 |

81 | 0.64 | 0.78 | 16.36 |

89 | 0.66 | 0.74 | 15.20 |

98 | 0.60 | 0.75 | 16.50 |

108 | 0.61 | 0.76 | 17.82 |

119 | 0.63 | 0.76 | 19.15 |

131 | 0.64 | 0.79 | 21.59 |

144 | 0.59 | 0.69 | 17.45 |

158 | 0.58 | 0.67 | 17.73 |

174 | 0.61 | 0.73 | 21.24 |

191 | 0.59 | 0.72 | 21.91 |

211 | 0.57 | 0.73 | 23.73 |

232 | 0.59 | 0.68 | 22.25 |

255 | 0.58 | 0.73 | 26.32 |

281 | 0.59 | 0.78 | 30.95 |

309 | 0.53 | 0.71 | 27.11 |

340 | 0.61 | 0.79 | 35.79 |

374 | 0.55 | 0.79 | 37.85 |

411 | 0.58 | 0.77 | 37.49 |

452 | 0.60 | 0.79 | 41.18 |

497 | 0.62 | 0.87 | 58.07 |

547 | 0.56 | 0.81 | 47.56 |

602 | 0.55 | 0.81 | 50.05 |

662 | 0.58 | 0.85 | 60.43 |

728 | 0.54 | 0.79 | 52.95 |

801 | 0.54 | 0.76 | 50.76 |

881 | 0.51 | 0.75 | 51.47 |

970 | 0.59 | 0.85 | 73.70 |

Building on MRA, we can add in lagged relationships to the analysis
using `mlra()`

and see not only whether there are scale-wise
variations in the *β* coefficients, but changes in these at
particular time lags.

```
<- logScale(scale.min = 10, scale.max = 1000, scale.ratio = 1.1)
scales <- mlra(x = sim1[,1], y = sim1[,2],order = 1, scales = scales, lags=100, direction='p')
mlra.out <- colorRampPalette(c("green", "red"))
colfunc ::matplot(mlra.out$betas, type='l', col = colfunc(49), ylab="Beta Coefficient", xlab='Lag', main="Multiscale Lagged Regression Analysis") graphics
```

The figure above shows that there is high *β* values across
scales only for lags near to 0. But, it’s difficult to see the
scale-wise variation in this type of plot. Another option is to use
`image()`

or `image.plot()`

to visualize the
results of the `mlra()`

function. This plot more clearly
shows a color gradient corresponding to the strength of the
/*b**e**t**a* coefficient across scales on the y-axis and
at increasing lags (x-axis).

```
<- 0:100
x <- scales
y image.plot(x, y, mlra.out$betas, axes=TRUE, legend.lab = "Beta Coefficient", ylab="Scale", xlab="Lag", main="Multiscale Lagged Regression Analysis")
```

`#contour(x, y, mlra.out,levels=seq(0,1,by=1),add=TRUE,col='black')`

So far we’ve demonstrated analyses on simulated data. Next we are
going to analyze a pair of empirical handmovement time series from a
dyad that are included in the `crqa`

package and discussed in
detail in Wallot et al. (2016). These data are hand-movement velocity
profiles.

```
require(crqa)
data(crqa)
head(handmovement)
```

```
## P1_TT_d P1_TT_n P2_TT_d P2_TT_n
## 1 0.22702423 0.26757616 0.006000000 0.000000000
## 2 0.47086091 0.23433310 0.008366600 0.006708204
## 3 0.14652304 0.16225289 0.006708204 0.006000000
## 4 0.40472830 0.02319483 0.009055385 0.010392305
## 5 0.07187489 0.04341659 0.007348469 0.006708204
## 6 0.10543719 0.06870953 0.005196152 0.007348469
```

```
plot(handmovement$P1_TT_d, type='l', main = "Dominant hand-movement velocity profiles of participants", xlab='time',ylab='movement velocity')
lines(handmovement$P2_TT_d, col=2)
```

Using the dominant hand variable for participant 1
`P1_tt_d`

and for participant two `P1_tt_d`

, we
are going to look at the fractal dynamics first. We are going to run
`mfdfa()`

and look at the q=2 (monofractal) scaling as well
as multifractal scaling.

```
# Participant 1
<- seq(-5,5, by=0.25)
q
<- logscale(scale_min = 16, scale_max = length(handmovement$P1_TT_d)/4, scale_ratio = 2)
scales
<- mfdfa(x = handmovement$P1_TT_d, q = q, order = 1, scales = scales)
mf.dfa.P1hand.out
# Participant 2
<- mfdfa(x = handmovement$P2_TT_d, q = q, order = 1, scales = scales)
mf.dfa.P2hand.out
# Examine the alpha exponent for q=2, which is equivalent to running DFA
$Hq[mf.dfa.P1hand.out$q==2] mf.dfa.P1hand.out
```

`## [1] 0.8788367`

`$Hq[mf.dfa.P2hand.out$q==2] mf.dfa.P2hand.out`

`## [1] 0.9174263`

For Participant 1, we observe a value of 0.88 and for Participant 2 a value of 0.92, which suggests both exhibit long-range correlations and the signals approximate pink noise. Next, want to examine the multi-fractal spectra.

```
# Create the plot
plot(mf.dfa.P1hand.out$h,mf.dfa.P1hand.out$Dh, type='b', lwd=1, lty=2, pch=19,ylim=c(-0.4,1),xlim=c(-.8,.8), ylab="D(h)", xlab="h", main= "Multifractal Spectrum for Dominant Hand Movements")
lines(mf.dfa.P2hand.out$h,mf.dfa.P2hand.out$Dh, type='b', pch=19,lwd=3, col="blue")
legend(-.85,1, legend = c("P1", "P2"), col=c("black", "blue"), lwd=3)
```

```
<- max(mf.dfa.P1hand.out$h) - min(mf.dfa.P1hand.out$h)
P1.W <- max(mf.dfa.P2hand.out$h) - min(mf.dfa.P2hand.out$h) P2.W
```

When comparing multi-fractal spectrum width (W),
*h*_{max} − *h*_{min},
we observe that both signals have a similar W. Specifically, Participant
1 W = NaN and Participant 2 W = NaN. From the figures and W estimates,
there does appear to be scale-wise variation in the scaling exponents.
However, a surrogate test could make this inference more robust.

Next, we are going to work with these hand movement time series from a dyad and try to examine the scale-wise cross-correlation between the time series. But first, let’s see if they are cross-correlated in general.

`ccf(handmovement$P1_TT_d, handmovement$P2_TT_d, lag.max = 500, main = "Cross-correlation function of P1 and P2 Hand Movements")`

It appears that there might be some lagged relationship between the
two signals. Now, we can run and examine the results of
`dcca()`

examining the scale-wise detrended cross-correlation
coefficients.

```
# Set scales
<- seq(15, 1000, by = 5)
scales # Note that a small amount of noise was added to these time series to avoid processor specific issues.
# This is not a necessary step for typical analyses
= handmovement$P1_TT_d + rnorm(1, 0,.001)
p1 = handmovement$P2_TT_d + rnorm(1, 0,.001)
p2 = dcca(x = p1, y = p2, order = 1, scales = scales)
dcca.out <- as.data.frame(dcca.out)
dcca.out # dcca.plot <- ggplot(data=dcca.out, aes(x=scales,y=rho)) + geom_point() +geom_line()
plot(dcca.out$scales, dcca.out$rho, type = 'b', pch = 16, xlab = 'Scale',
ylab = expression(rho))
```

`# dcca.plot`

At smaller scales, *ρ* approximates zero. However, at
increasing scale sizes the trend of the *ρ* estimates become
negative exhibit quite large fluctuations.

Building on the above DCCA results, we use `mra()`

to
determine if can we use the the scale-wise fluctuations in Particiapnt
2’s hand movements to predict those in Participant 1.

```
<- seq(15, 1000, by = 5)
scales = handmovement$P1_TT_d + rnorm(1, 0, .001)
p1 = handmovement$P2_TT_d + rnorm(1, 0, .001)
p2 <- as.data.frame(mra(x = p1, y = p2, order = 2, scales = scales))
mra.out <- ggplot(data=mra.out, aes(x=scales,y=betas)) + geom_point() +geom_line()
mra.plot mra.plot
```

In the table below, we filter out those *β* coefficients,
whose corresponding t-values are greater than +/- 1.96 to provide an
index of how many scales there are where Participant 2’s hand movements
are significant predictors for Participant 1’s hand movements.

```
<- subset(mra.out, abs(mra.out$t_observed) > 1.96)
mra.out.sig ::kable(mra.out.sig, format='html', digits =2,align='c',caption = "Output from MRA on Handmovement Data") knitr
```

scales | betas | r2 | t_observed | |
---|---|---|---|---|

38 | 200 | -0.12 | 0.95 | -7.97 |

39 | 205 | -0.12 | 0.85 | -4.37 |

42 | 220 | -0.12 | 0.69 | -2.79 |

43 | 225 | -0.10 | 0.70 | -2.82 |

44 | 230 | -0.17 | 0.67 | -4.03 |

45 | 235 | -0.11 | 0.67 | -2.66 |

46 | 240 | -0.10 | 0.74 | -3.57 |

47 | 245 | -0.08 | 0.70 | -2.25 |

48 | 250 | -0.20 | 0.74 | -6.76 |

49 | 255 | -0.09 | 0.75 | -3.75 |

50 | 260 | -0.15 | 0.64 | -4.54 |

52 | 270 | -0.21 | 0.54 | -4.72 |

55 | 285 | -0.08 | 0.60 | -2.46 |

57 | 295 | -0.14 | 0.42 | -3.32 |

59 | 305 | 0.07 | 0.52 | 2.03 |

61 | 315 | -0.15 | 0.44 | -3.89 |

62 | 320 | -0.09 | 0.44 | -2.55 |

63 | 325 | -0.15 | 0.39 | -3.40 |

64 | 330 | -0.16 | 0.43 | -4.03 |

65 | 335 | -0.14 | 0.47 | -3.95 |

66 | 340 | -0.09 | 0.46 | -2.82 |

67 | 345 | -0.09 | 0.38 | -2.51 |

68 | 350 | -0.18 | 0.45 | -6.07 |

69 | 355 | -0.15 | 0.44 | -4.55 |

72 | 370 | -0.25 | 0.35 | -5.38 |

73 | 375 | -0.08 | 0.42 | -2.26 |

82 | 420 | -0.09 | 0.24 | -2.48 |

83 | 425 | -0.15 | 0.31 | -4.33 |

84 | 430 | -0.08 | 0.33 | -2.30 |

88 | 450 | -0.12 | 0.17 | -2.39 |

89 | 455 | -0.25 | 0.21 | -5.08 |

90 | 460 | -0.23 | 0.22 | -5.36 |

91 | 465 | -0.16 | 0.22 | -4.55 |

92 | 470 | -0.11 | 0.21 | -3.31 |

93 | 475 | -0.09 | 0.19 | -2.40 |

94 | 480 | -0.14 | 0.19 | -3.16 |

95 | 485 | -0.32 | 0.22 | -6.37 |

96 | 490 | -0.41 | 0.29 | -9.29 |

97 | 495 | -0.38 | 0.30 | -9.45 |

98 | 500 | -0.26 | 0.23 | -6.10 |

99 | 505 | -0.14 | 0.18 | -3.18 |

100 | 510 | -0.14 | 0.19 | -3.24 |

101 | 515 | -0.14 | 0.18 | -3.36 |

102 | 520 | -0.15 | 0.18 | -4.03 |

103 | 525 | -0.19 | 0.19 | -5.55 |

104 | 530 | -0.25 | 0.21 | -7.50 |

105 | 535 | -0.29 | 0.24 | -8.88 |

106 | 540 | -0.31 | 0.23 | -9.09 |

107 | 545 | -0.36 | 0.23 | -9.43 |

108 | 550 | -0.34 | 0.20 | -8.47 |

109 | 555 | -0.30 | 0.16 | -6.88 |

110 | 560 | -0.28 | 0.13 | -6.01 |

111 | 565 | -0.32 | 0.13 | -6.18 |

112 | 570 | -0.34 | 0.13 | -6.42 |

113 | 575 | -0.37 | 0.14 | -6.94 |

114 | 580 | -0.47 | 0.17 | -8.74 |

115 | 585 | -0.58 | 0.24 | -11.28 |

116 | 590 | -0.60 | 0.27 | -12.64 |

117 | 595 | -0.54 | 0.27 | -12.50 |

118 | 600 | -0.52 | 0.28 | -13.11 |

119 | 605 | -0.46 | 0.26 | -12.24 |

120 | 610 | -0.43 | 0.25 | -12.08 |

121 | 615 | -0.37 | 0.23 | -11.06 |

122 | 620 | -0.32 | 0.19 | -9.38 |

123 | 625 | -0.33 | 0.18 | -9.02 |

124 | 630 | -0.40 | 0.19 | -9.73 |

125 | 635 | -0.55 | 0.25 | -12.37 |

126 | 640 | -0.63 | 0.29 | -14.14 |

127 | 645 | -0.69 | 0.33 | -15.91 |

128 | 650 | -0.68 | 0.33 | -16.31 |

129 | 655 | -0.62 | 0.31 | -15.37 |

130 | 660 | -0.55 | 0.27 | -13.58 |

131 | 665 | -0.44 | 0.19 | -10.13 |

132 | 670 | -0.31 | 0.12 | -6.66 |

133 | 675 | -0.23 | 0.08 | -4.60 |

134 | 680 | -0.19 | 0.07 | -3.74 |

135 | 685 | -0.18 | 0.07 | -3.57 |

136 | 690 | -0.19 | 0.07 | -3.67 |

137 | 695 | -0.21 | 0.07 | -4.22 |

138 | 700 | -0.25 | 0.09 | -5.38 |

139 | 705 | -0.27 | 0.10 | -6.34 |

140 | 710 | -0.29 | 0.13 | -7.44 |

141 | 715 | -0.31 | 0.15 | -8.48 |

142 | 720 | -0.28 | 0.15 | -8.35 |

143 | 725 | -0.26 | 0.14 | -7.94 |

144 | 730 | -0.23 | 0.12 | -7.11 |

145 | 735 | -0.20 | 0.11 | -6.21 |

146 | 740 | -0.17 | 0.10 | -5.27 |

147 | 745 | -0.16 | 0.09 | -4.56 |

148 | 750 | -0.18 | 0.09 | -4.83 |

149 | 755 | -0.23 | 0.10 | -5.73 |

150 | 760 | -0.27 | 0.11 | -6.97 |

151 | 765 | -0.30 | 0.13 | -7.93 |

152 | 770 | -0.29 | 0.13 | -7.95 |

153 | 775 | -0.29 | 0.13 | -8.02 |

154 | 780 | -0.32 | 0.13 | -8.56 |

155 | 785 | -0.35 | 0.14 | -9.03 |

156 | 790 | -0.39 | 0.14 | -9.33 |

157 | 795 | -0.40 | 0.14 | -9.37 |

158 | 800 | -0.39 | 0.13 | -8.82 |

159 | 805 | -0.35 | 0.11 | -7.83 |

160 | 810 | -0.36 | 0.11 | -7.78 |

161 | 815 | -0.39 | 0.12 | -8.34 |

162 | 820 | -0.42 | 0.13 | -9.17 |

163 | 825 | -0.46 | 0.14 | -9.99 |

164 | 830 | -0.50 | 0.16 | -11.13 |

165 | 835 | -0.52 | 0.18 | -12.12 |

166 | 840 | -0.53 | 0.19 | -12.82 |

167 | 845 | -0.52 | 0.19 | -13.03 |

168 | 850 | -0.50 | 0.19 | -12.99 |

169 | 855 | -0.48 | 0.19 | -12.95 |

170 | 860 | -0.45 | 0.18 | -12.62 |

171 | 865 | -0.41 | 0.16 | -11.85 |

172 | 870 | -0.37 | 0.15 | -10.95 |

173 | 875 | -0.34 | 0.13 | -10.09 |

174 | 880 | -0.31 | 0.12 | -9.41 |

175 | 885 | -0.29 | 0.11 | -8.86 |

176 | 890 | -0.27 | 0.10 | -8.22 |

177 | 895 | -0.26 | 0.09 | -7.69 |

178 | 900 | -0.26 | 0.09 | -7.48 |

179 | 905 | -0.27 | 0.09 | -7.75 |

180 | 910 | -0.30 | 0.10 | -8.29 |

181 | 915 | -0.33 | 0.11 | -9.06 |

182 | 920 | -0.36 | 0.13 | -10.15 |

183 | 925 | -0.42 | 0.16 | -11.84 |

184 | 930 | -0.50 | 0.20 | -14.20 |

185 | 935 | -0.58 | 0.26 | -17.01 |

186 | 940 | -0.66 | 0.32 | -20.14 |

187 | 945 | -0.71 | 0.38 | -23.32 |

188 | 950 | -0.74 | 0.44 | -26.35 |

189 | 955 | -0.74 | 0.47 | -28.26 |

190 | 960 | -0.71 | 0.48 | -28.60 |

191 | 965 | -0.67 | 0.45 | -26.76 |

192 | 970 | -0.63 | 0.42 | -25.38 |

193 | 975 | -0.60 | 0.39 | -23.74 |

194 | 980 | -0.57 | 0.35 | -21.71 |

195 | 985 | -0.55 | 0.30 | -19.49 |

196 | 990 | -0.52 | 0.26 | -17.27 |

197 | 995 | -0.46 | 0.20 | -14.50 |

198 | 1000 | -0.44 | 0.17 | -12.80 |

Lastly, let’s take a look at the MLRA plot. Below we can see that most of the highest beta coefficients are at very short lags; however, there is some considerable variability around the scales.

```
<- logScale(scale.min = 10, scale.max = 1000, scale.ratio = 1.1)
scales <- mlra(x = handmovement$P1_TT_d, y = handmovement$P2_TT_d,order = 1, scales = scales, lags=100, direction='p')
mlra.out.emp <- 0:100
x <- scales
y image.plot(x, y, mlra.out.emp$betas, axes=TRUE, legend.lab = "Beta Coefficient", ylab="Scale", xlab="Lag", main="Multiscale Lagged Regression Analysis Hand Movements")
```

`#contour(x, y, mlra.out,levels=seq(0,1,by=1),add=TRUE,col='black')`

- Ihlen, E. A. F. (2012). Introduction to Multifractal Detrended Fluctuation Analysis in Matlab. Frontiers in Physiology, 3. https://doi.org/10.3389/fphys.2012.00141
- Kantelhardt, J. W., Koscielny-Bunde, E., Rego, H. H. A., Havlin, S., & Bunde, A. (2001). Detecting long-range correlations with detrended fluctuation analysis. Physica A: Statistical Mechanics and Its Applications, 295(3), 441–454. https://doi.org/10.1016/S0378-4371(01)00144-3
- Kantelhardt, J. W., Zschiegner, S. A., Koscielny-Bunde, E., Havlin, S., Bunde, A., & Stanley, H. E. (2002). Multifractal detrended fluctuation analysis of nonstationary time series. Physica A: Statistical Mechanics and Its Applications, 316(1), 87–114. https://doi.org/10.1016/S0378-4371(02)01383-3
- Kristoufek, L. (2013). Mixed-correlated ARFIMA processes for power-law cross-correlations. Physica A: Statistical Mechanics and Its Applications, 392(24), 6484–6493. https://doi.org/10.1016/j.physa.2013.08.041
- Peng, C.-K., Buldyrev, S. V., Havlin, S., Simons, M., Stanley, H. E., & Goldberger, A. L. (1994). Mosaic organization of DNA nucleotides. Physical Review E, 49(2), 1685–1689. https://doi.org/10.1103/PhysRevE.49.1685
- Podobnik, B., & Stanley, H. E. (2008). Detrended Cross-Correlation Analysis: A New Method for Analyzing Two Nonstationary Time Series. Physical Review Letters, 100(8), 084102. https://doi.org/10.1103/PhysRevLett.100.084102
- Wallot, S., Mitkidis, P., McGraw, J. J. and Roepstorff, A. (2016). Beyond synchrony: joint action in a complex production task reveals beneficial effects of decreased interpersonal synchrony. PloS one, 11(12), e0168306.