Assessing Time Series Correlation Significance: A Parametric Approach with Application to Physiological Signals

Correlation coefficients play a pivotal role in quantifying linear relationships between random variables. Yet, their application to time series data is very challenging due to temporal dependencies. This paper introduces a novel approach to estimate the statistical significance of correlation coefficients in time series data, addressing the limitations of traditional methods based on the concept of effective degrees of freedom (or effective sample size, ESS). These effective degrees of freedom represent the independent sample size that would yield comparable test statistics under the assumption of no temporal correlation. We propose to assume a parametric Gaussian form for the autocorrelation function. We show that this assumption, motivated by a Laplace approximation, enables a simple estimator of the ESS that depends only on the temporal derivatives of the time series. Through numerical experiments, we show that the proposed approach yields accurate statistics while significantly reducing computational overhead. In addition, we evaluate the adequacy of our approach on real physiological signals, for assessing the connectivity measures in electrophysiology and detecting correlated arm movements in motion capture data. Our methodology provides a simple tool for researchers working with time series data, enabling robust hypothesis testing in the presence of temporal dependencies.


Introduction
Correlation coefficients, such as Pearson's or Spearman's, are fundamental tools for assessing linear relationships between random variables.Although originally not designed for dependent samples, these coefficients are widely used with time series in diverse fields.Notable examples can be found in neuroimaging, where Pearson's correlation coefficient is used to construct functional brain networks from functional Magnetic Resonance Imaging (fMRI) data (Fox et al., 2005;Van Dijk et al., 2010), electrophysiology data (Zhang et al., 2021;Naira et al., 2019;Ji et al., 2019;Zhong et al., 2022) or to relate brain signals to other behavioral parameters, such as movement data (Lu et al., 2021).In this paper, we address the estimation of the statistical significance of correlation coefficients for time series.
Fisher's variance-stabilizing transformation is useful to obtain simple bounds and significance of a correlation coefficient r.Under the null hypothesis that two sets of n independent data points are uncorrelated, the null distribution of the Fisher-transformed variable z = arctanh(r) is approximately normally distributed with mean 0 and standard deviation 1/ √ n − 3.In other words, we have In real time-series data, the assumption of pairwise independence underlying traditional correlation coefficient tests often falters due to dependencies between consecutive observations.
When data exhibit temporal dependencies, Eq. ( 1) is overconfident and fails to account for the loss of degrees of freedom due to temporal correlations.A pivotal solution, pioneered by Bartlett in 1935(Bartlett, 1935), adjusts for this change of degrees of freedom by introducing a number of "effective degrees of freedom", or "effective sample size" (ESS), representing the independent sample size producing comparable test statistics.Under the null hypothesis that two sets of n dependent data points, with ESS ν (ν ≤ n), are uncorrelated, we can simply rewrite Eq. ( 1) to obtain the approximate distribution as: This enables hypothesis testing even when actual samples show temporal dependencies.
The approach pioneered by Bartlett, that had various extensions, e.g., (Bartlett, 1946;Bayley and Hammersley, 1946;Quenouille, 1947;Pyper and Peterman, 1998;Afyouni et al., 2019), relies on estimating the ESS from the sum of the product of the autocorrelation functions (ACFs) of the time series.Let ρ k and γ k be the autocorrelation at lag k of two time series, resp.x and y, of n samples each.The ESS of the correlation coefficient between x and y is (Quenouille, 1947): While this estimator of the ESS has been widely adopted, it requires computing the sample ACF which can be computationally demanding and can also yield inaccurate estimates of the ESS due to accumulating noise (in the sum of Eq. ( 3)).
In this work, we assume a parametric form for the ACF to simplify the computation and estimation of the ESS.In particular, we show that the sum in Eq. (3) converges to an integral which can be analytically solved under a Gaussian (Laplace) approximation.The resulting expression relates the ESS to the average second spectral moment, also called roughness, of the pair of time series to correlate.Importantly, the roughness of a series can be estimated from the variance of its temporal derivatives.This scaffolds the central result of this work: for two processes x and y of length n with temporal derivatives ẋ and ẏ, the number of effective degrees of freedoms ν is approximately: We show that this formula yields accurate estimation of the ESS and effectively provides a computationally effective method to evaluate the significance of correlation coefficients.
This article is organized as follows.First, we present the key contributions of our work, which encompass the derivation of the asymptotic expression for ESS, a straightforward approximation for Gaussian autocorrelation, and the introduction of an estimator based on the variance of the temporal derivatives of the processes.Then, we validate our approach through numerical experiments, demonstrating its effectiveness and its computational efficacy.Finally, we apply our methodology to real data to (i) highlight the significance of power-based connectivity measures in electrophysiology, and (ii) distinguish correlated and uncorrelated arm captured motion data.

A parametric estimator of the effective sample size
Our work builds on estimating the ESS using a Gaussian approximation to the ACF.In this section we show that the ESS converges to an integral; then that a Laplace approximation of this integral relates the ESS to the signals roughness; and finally, that common roughness estimators can be used to build an estimator of the ESS.
This work focuses on estimating the ESS for the correlation of smooth, wide-sense stationary (WSS) stochastic processes.A smooth WSS stochastic process x is defined by convolving white noise w with a smooth (C ∞ ), square-integrable function K : R → R + , in other word: For clarity and brevity, we reuse the same notation throughout this paper.By default x and y are two smooth WSS stochastic processes with zero mean, variances resp.σ x and σ y , and autocorrelations functions ρ and γ.For simplicity, we restrict our result section to the case where both x and y have the same autocorrelation, i.e., ρ = γ.This restriction is moreover motivated in Section 2.2.

Asymptotic expression for the ESS
This section introduces an asymptotic form for the ESS.We consider Eq. ( 3) for the ESS in the case of infinitely large n and infinitely small sampling interval of the ACF.Under these conditions, the sum in Eq. ( 3) converges to This yields a new asymptotic ESS expression, ν ∞ : The resulting integral formulation is pivotal in deriving analytical expressions for the ESS of stochastic processes with known ACF.As for processes with known autocorrelation function, we can directly evaluate Eq. ( 7).In the next section, we develop a generic approximation based on a Laplace approximation of the integral.

Closed-form expression under a Laplace approximation
A general approach to approximate the integral in the asymptotic ESS expression is to use Laplace's method.In this section, we derive a general approximation to asymptotic ESS (Eq.( 7)), based on a Laplace approximation of the integral (Eq.( 6)).In brief, Laplace's method is used to approximate the integral of a arbitrary function with a unique global maximum using a simpler, Gaussian function around its peak (Penny et al., 2011).Typically in physiological signals, ACFs have a mode at lag 0 and decrease to 0 for large lags, although they might exhibit long-range correlations.Thus, the area under the curve of the product of two such ACFs is expected to be contained around 0, with an attenuation of long-range correlations.This justifies using Laplace's method to approximate the asymptotic ESS expression.
To derive the expression under Laplace approximation, note that the second-order expansion of the ACFs ρ and γ around their mode at τ = 0 is similar to processes with Gaussian ACFs where ρ ′′ (0) and γ ′′ (0) are the second spectral moment of the processes, i.e., the second-order derivatives of the ACFs at 0. The second spectral moment is a universal measure of "roughness" in the literature of stochastic processes (Friston et al., 2008;Cox and Miller, 2017).From this Gaussian form, we can explicitly evaluate the integral in Eq. ( 7): By substituting Eq. ( 9) in Eq. ( 7), the ESS under a Laplace approximation is To summarise, we can use a Laplace (Gaussian) approximation to approximate the ESS.This is equivalent to fitting the ACF using a one-parameter family of Gaussian functions parameterised by their second spectral moment, or roughness.This allows to derive a simple expression for the ESS that only involves the second spectral moment of each process.

Estimating the second spectral moment
A well-established result from stochastic process theory is that the roughness can be conveniently estimated from the variance of the first-order temporal derivatives of the process (Cox and Miller, 2017;Adler, 2010;Worsley, 1996).In this section, we leverage this result to construct an estimator of the asymptotic ESS.
Let x and y be two processes with autocorrelation functions ρ and γ and temporal derivatives  10), we obtain an estimator of the ESS: By construction, this new parametric estimator for the ESS depends only on the variance of the temporal derivatives of the time series.As such, its consistency and unbiasedness derive directly from those of the variance estimators (from direct application of the continuous mapping theorem).
Remark 2. The second spectral moment is also related to the expected number of zero-crossings of a process by Rice's formula (Rice, 1944).For didactic purposes, we show an alternative ESS estimator based on Rice's formula in Appendix A.

Numerical validation
We conduct numerical experiments on Gaussian processes with Gaussian autocorrelation functions (GPGA) to assess the validity and limitations of our approach.The numerical results are organised in four parts.First, we evaluate the quality of fitting the squared autocorrelations from the process roughness.Second, we compare the quality of the estimates of the ESS, corresponding to the integral of the squared ACFs, with existing methods.Then, we compare the statistics yielded by our method with existing ones.Finally, we investigate the computational performance of our approach.

Assessment of parametric autocorrelation function estimates
The method we present leverages the analytical evaluation of the ESS of a pair of stochastic processes from their roughness, estimated from the variance of the temporal derivatives.Here, we verify whether this method can accurately estimate the ACF of the process, when the ACF is known to be Gaussian.
We generate 2000 sample paths from GPGA with different roughness levels.We then retrieve the squared autocorrelation using the variance of the derivatives.To sample a sample path from a GPGA with roughness r, we convolve a unit white Gaussian noise sample with a Gaussian kernel of variance 2/r as per (Friston et al., 2008).We then apply Eq. ( 4) to determine its second spectral moment from the variance of the process temporal derivatives.
Figure 1 illustrates that low roughness results in significant biases in squared sample autocorrelation functions, leading to ESS underestimation.This bias becomes clear in Figure 1 where the grey curves diverge significantly from their expected near-zero values.This figure underscores the merit of predetermining the functional form of the ACF, thereby filtering out spurious long-term correlations.Furthermore, our method consistently aligns more closely with the theoretical ACF relatively to the direct sample squared ACF.
In conclusion, explicitly formulating the functional form of the ACF filters out errant longterm correlations otherwise amplified by traditional ESS estimators based on sample autocorrelations.This leads to more accurate and consistent ACF estimates across varying roughness scales.However, we note that the variance of the squared ACF appears to be inversely proportional to roughness.This issue is further investigated in subsequent sections.Figure 1: Sample paths and squared autocorrelation functions for roughness values ρ ′′ (0).Rows one and two display, respectively, samples of Gaussian processes for different roughness values and their squared autocorrelation functions.Notably, the squared ACF plays a crucial role in ESS estimation.The plots showcase various ACF estimations: theoretical (plain yellow), sample (plain grey), average of sample (dotted black), Gaussian fit based on estimated roughness (plain red), and Gaussian fit based on average roughness (dashed blue).The bias in sample ACFs, evident in the grey curves, increases as roughness decreases due to increasing random long-range correlations.On the other hand, the Gaussian fits (plain red) remain unbiased at large lags, ensuring the ACF retains its Gaussian form.

Assessment of parametric ESS estimates 1
This section deals with the influence of the roughness estimator's variance on the variance of 2 the estimated sum of squared autocorrelation and that of the ESS.We also examine how these 3 relate to the length of series.

4
We sampled 1000 sample paths from GPGA, adjusting roughness between 10 −6 and 1.Each 5 sample path's roughness is then estimated.Using this roughness estimate, we compute the 6 ESS using Eq. ( 11).Furthermore, we determine the ACF of the process using the inverse 7 Fourier transform of the sample path's power spectral density (PSD).The latter is estimated 8 either through the FFT or the Welch periodogram with a 256-point window having 128-point 9 overlap (Welch, 1967).The ESS is then obtained from the ACF using Eq.(3).We replicate this 10 procedure for series lengths of 500, 1000, and 2000 points.For comparison, the ESS is normalized From Figure 2(a), it is clear that roughness estimates share a similar trend.For minuscule 1 roughness values, bias emerges.Extended series better approximate low roughness, hinting that 2 biases result from series that are too short relative to their roughness.As expected from their 3 relationship (Eq.( 7)), ESS factors derived from roughness mirror this trend.In summary, our numerical analyses confirm the paramountcy of a parametric autocorrelation function, where the unique parameter is estimated via the process's second spectral moment.
This method circumvents biases introduced by random long-range dependencies, resulting in more accurate sample p-values for correlation coefficients.In subsequent sections, we advocate for GPGA's aptness in real-world applications, such as neuroimaging data connectivity analysis and movement trajectory correlations.

Comparison of computational performance
In this section, we evaluate the performance gain of our approach.The main difference with existing methods is that we avoid explicit evaluation of the ACF by using a parametric ACF, which only requires to estimate the average series roughness.
To evaluate the impact of series length on the computation time and speedup, we generate pairs of sample paths from a GPGA process with roughness 10 −3 having varying length, from 100 to 10 7 points with 10 fold increments.For each pairs of paths, we compute the ESS using FFT-based and Welch-based computation of the ACFs, and with our approach.The FFT-based ACF is obtained using numpy's fft and ifft functions (Harris et al., 2020).The Welch-based ACF is computed using scipy's welch function from the signal package with a window length of 256 points (100 points for the 100 points series) and numpy's ifft (Virtanen et al., 2020).
Our approach uses numpy's diff and var to compute the variance of the temporal derivatives of the process.The results are obtained on a laptop equipped with an 12-core Intel Xeon W-10855M CPU and 32Gb of RAM running Ubuntu 22.04 and Python 3.8.For each approach, we computed the timings from 100 loops and present the mean and standard deviation across 7 different runs.
We report the results on Speedup of our approach, obtain as the average computation time of each method divided by the average computation time of the proposed method.For instance, a speedup of 5 means that our approach divides by 5 the average computation time.
tational advantage of the proposed approach.In average across different number of points, the 1 proposed approach is 7.4 times faster than an FFT-based approach.This speedup is negligible 2 for short time series (resp.1.3 and 1.6 for 100 and 1000 points) but significant for long time se-3 ries (10.2 to 13.4 for 10 5 to 10 7 ).The speedup from Welch-based approach is relatively constant 4 across time series length due to the fixed window length, with an average speedup of 5.1.

5
To sum-up, our proposed approach, based on the temporal derivatives of the time series, gives a significant speedup as compared to other approaches based on the computation of the 7 ACF.This computational argument motivates using a parametric Gaussian form of the ACF in 8 applications where timing is critical or where large amounts of data have to be processed; for 9 instance, when evaluating correlations of brain activity between a large number of brain regions.Step 4.

Test against the null hypothesis
Step 1.

Compute the time-frequency decomposition
Step 2.
Extract the 4th roots of power time series p-value Step 3. uate the functional connectivity between brain regions.A standard approach, power-based connectivity analysis, measures functional connectivity by evaluating the correlation between the time-frequency power of the signal at two different locations (Cohen, 2014).Several methods exist to perform power-based connectivity analysis, but the correlation between wavelet-based time-frequency representations is particularly interesting within the scope of this work.

Compute the correlation coefficient and the ESS
The time course of power can be obtained using a continuous wavelet transform with a Morlet wavelet.A Morlet wavelet combines a Gaussian kernel with a complex sinusoid at a specified frequency.Therefore, a continuous wavelet transform with a Morlet wavelet introduces a form of temporal smoothing with a Gaussian kernel.Intuitively, if the width of the Gaussian smoothing kernel is large enough relative to the autocorrelation of the signal, the resulting time series will have roughly a Gaussian autocorrelation.In this case, one can consider testing the strength of the correlation between two power time series under the null hypothesis that the signals are two uncorrelated GPGA.The overall procedure is illustrated in Figure 5.This approach allows a more efficient analysis of connectivity significance without relying on computationally demanding methods.
This section aims at gaining some insight on whether the null hypothesis of uncorrelated GPGA processes might yield appropriate statistics.In contrast to the numerical experiments section, we cannot sample from the null distribution of correlation coefficient between bandpower time series.Thus, we can only investigate a few characteristics of the data to assess the validity of our approach.Indeed, this section echoes the first step of statistical analysis in real world application, where statisticians have no access to the data generative distribution and shall assess the suitability of their model.

Data presentation and preprocessing
We analyse a dataset containing EEG recordings of a subject listening to continuous, naturalistic speech (Di Liberto et al., 2015;Crosse et al., 2016).The subject listened to a classic work read by an English speaker.EEG signals were recorded from a single participant with a Biosemi system having 128 channels and a sampling rate of 512 Hz.Prior to analysis, the data has been filtered between 1 and 15 Hz and downsampled to 128Hz.We then computed the time-frequency representation of the EEG signals using a continuous wavelet transform.We used Morlet wavelet with 7 cycles to produce time-frequency representation with a sufficient temporal smoothness to fall under our Gaussian assumptions.Finally, we take the fourth root of the power to make the data marginally Gaussian (Hawkins and Wixley, 1986).
Remark 3. Note that selecting the number of cycles in the wavelet effectively selects the degree of smoothness of the signal.This shows that the number of cycles is related to the ESS.In Appendix B, we give a formula giving the ESS as a function of the number of cycles.

Adequacy of GPGA assumptions for wavelet-based EEG power
Here, we show that our approach yieds appropriate statistics for power-based connectivity analysis.We first look at whether the signals have an approximately Gaussian autocorrelation and whether our approach provides a good fit of the squared autocorrelation function.We compare results given by FFT-based and Welch-based approaches with our approach, and analyse the results in the light of our numerical results.Finally, we compare the ESS and 97.5% quantile given by all three approaches with real data and sample GPGA sample paths with matching roughness.
Figure 6 summarizes the properties of the data.In particular, we show a time-frequency plot (Fig. 6 (a) and (b)) and some samples power time series at different frequencies.Intuitively, the time course of power at lower frequencies is smoother than for higher frequencies.This is confirmed by looking at the roughness of the signals across channels, whose mean increase with frequency (Fig. 6 (d)).In addition, we display (Fig. 6 (c)) the marginal distribution of the signals  across channels to inspect its normality -a key element in mandating using Pearson's correlation 1 coefficient instead of Spearman's rank correlation.Finally, we inspect (Fig. 6 (e), (f), and (g)) the squared ACF at 4 Hz, 10 Hz, and 14 Hz, and compare it against the Gaussian fit obtained 3 from the signals roughness.We observe that the squared ACF is roughly Gaussian, despite a 4 slower decrease (lower kurtosis) which might come from autocorrelations inherent to the nature 5 of the signal.Overall, these results suggest that GPGA might provide a good approximation 6 for the data and that our approach gives appropriate statistics.

Assessment of ESS and quantiles for wavelet-based EEG power
To confirm these results, we evaluate the ESS factor for all possible between-channel pairs connectivity at each given frequency.In addition, for each signal we generate a random GPGA sample path with matching length and roughness, and compute the ESS factor for all pairs of random sample paths.This gives us a hint on how the ESS factor would behave if the sample paths where effectively sampled from the null hypothesis, i.e., GPGA with matching length and roughness.We complement this analysis by looking at the 97.5% quantile of both real and simulated data, under the null hypothesis.Any correlation coefficient greater than the 97.5% quantile -or lower than its opposite -would be effectively considered as significant with p < 0.05.The 97.5% quantile is computed by the inverse Fisher transform of the 97.5% quantile of the Gaussian distribution: The ESS factor and 97.5% quantile for both real and simulated data are presented on Figure 7.
We observe that the ESS factor for real EEG signal increases with frequency.This is expected as the roughness of power time series increases with the frequency, as shown in Figure 6.This increased ESS factor results in a decreased 97.5% quantile.Indeed, as the effective sample size increases and statistical power increases, we can arbitrate on the rejection of the null hypothesis at lower values of correlation coefficient.This is conform to the role played by the ESS in Eq. ( 12).
Interestingly, we see that both FFT-based and Welch-based approaches give lower values of ESS, and higher 97.5% quantile value, than our proposed approach.This effect is the strongest at 10 Hz and 12 Hz.Comparing each figure with its analogous generated under GPGA assumptions, we see that under GPGA assumptions all three methods yield approximately the same values of ESS factor and 97.5% quantile.Because these differences are stronger at 10 Hz and 12 Hz, i.e., in the well-identified alpha band, we hypothesis that task-related modulations of the power task cause variations in the autocorrelation of the power; which induce variability and shift in the ESS factor and related quantiles.Note that the size of this effect on the 97.5% quantile is relatively small (less than 0.02).We argue that this is negligible in most cases, especially given the variability of the sample estimates of the roughness or ACFs.Gaussian.Thus, we aim at evaluating our proposed approach on less 'well-behaved' biological approach retrieves the statistics of correlation coefficients under the null hypothesis that signals are GPGA.For this particular type of process, our approach outperforms classical approaches, yielding robust estimate of the statistics under a large range of roughness.In addition, our approach shows relatively higher computational performances compared to alternative methods.
In a second part, we show that some data issued from biological signals satisfy GPGA assumptions.Our approach seem adequate to test for significant correlations between brain regions, when used with power-based connectivity analysis of electrophysiological data.Our method also applies well to random joint trajectories, as measured from a motion capture system.
Overall, obtained results suggest that our approach yields accurate and robust statistics under assumptions that can be found in real biological signal.More generally, we claim that using a parametric form of autocorrelation can yield more accurate statistics, without being overly restrictive on the form of the process.In particular, we think that more advanced methods for parameter estimation could allow better estimation of the autocorrelation function.However, as often with parametric approaches, there is no one-size-fits-all form of the autocorrelation function, and the statistician holds the responsibility to evaluate the adequacy of the model for the data under study.
In addition to correlation coefficients, our approach can be straightforwardly extended to linear regression.A pivotal assumption of linear regression is the independence of residuals.
An assumption often violated in time series data due to inherent temporal dependencies.The parametric ESS estimators derived here can be seamlessly adapted to assess the significance of regression coefficients in linear regression models applied to time series data.Thus, our approach not only fortifies the foundation of correlation analysis in time series but also extends its robustness to the broader landscape of linear regression.
of cycles in Morlet wavelets can be adjusted, providing a direct assessment of the ESS for the correlations.If two series, derived from the same number of cycles but different frequencies f 1 and f 2 , are correlated, their ESS is: For f 1 = f 2 = f , the ESS is simply ν ∞ = √ πnf /N c .Importantly, the chosen number of cycles determines the ESS of power correlations.

Figure 2
Figure 2: a) Estimated roughness as a function of the process roughness, for series of 500 points (dark blue), 1000 points (light blue), and 2000 points (yellow).b-d) Estimated ESS factor as a function of the process true ESS, computed from Eq. (10), for series of 500 points (dark blue), 1000 points (light blue), and 2000 points (yellow) for the FFT-based (b),Welch-based (c), and proposed approaches (d).Shaded areas correspond to the 95% confidence intervals.

4Figure 3 :6
Figure3: Estimated probability against empirical probability for the different approaches.The first row illustrates a roughness of 10 −2 , while the second showcases 10 −4 .Each column represents different estimation methods: FFT-based (green),Welch-based (red), and our method (blue).

Figure 4 .Figure 4
Figure4: a) Comparison of the computation time for varying number of points from 100 to .b) Speedup of our approach, obtain as the average computation time of each method divided by the average computation time of the proposed method.For instance, a speedup of 5 means that our approach divides by 5 the average computation time.
) captures the electric potentials on the scalp.Its low-cost and 15 high-temporal-resolution make EEG a particularly appealing neuroimaging technique to eval-

FOIFOIFigure 5 :
Figure 5: Procedure to assess the connectivity of two EEG channels at a particular frequency of interest (FOI).

Figure 6 :
Figure 6: (a) Time-frequency plot of the fourth root of power for channel A9.(b) Sample fourthroot of power time series at different frequencies.(c) Marginal distribution of fourth-root of power time series across channels.(d) Roughness across channels at different frequency.(e-g) Squared ACF at 4 Hz, 10 Hz, and 14 Hz.Gray lines indicate different channels, the black line indicates the mean squared ACF across channels, the dotted blue line indicates the mean fit using the proposed approach, and the filled blue area shows the area that is integrated in the denominator of the ESS expression. 7

Figure 7 :2
Figure 7: (a) Distribution of the ESS factor evaluated with the FFT-based (green), the Welch-based (red), and our approach (blue) for pair of power time series at a given frequency.(b) Similar to (a) for random GPGA sample paths with roughness matching that of every pair of signals at a given frequency.(c) Distribution of the 97.5% quantile obtained from the ESS of each pair of real time series.(d) Similar to (c) for random GPGA sample paths with matching roughness.