The analysis of resting-state surface EEG data, collected in a single session from each participant under EO and EC states, was conducted separately for each condition. The initial 10 s of each condition’s recording were discarded, and the subsequent 4-min segment was windowed and divided into non-overlapping 6-s segments for filtering. We determined this epoch length based on literature findings: In the most recent functional connectivity analyses, to estimate graph theoretic connectivity indices from simulated resting state 128-channel EEG data, mean squared COH, iCOH and wPLI were used for epoch lengths of 2 s, 4 s, 6 s, and 8 s and the number of epochs was varied to be 20, 40, max. As a result of detailed comparisons, it was reported that the most suitable connectivity estimates were obtained for 40 epochs of 6 s length (Miljevic et al. 2025). In previous studies calculating connectivity indices from eyes-closed resting-state 256-channel EEG data collected from the same participants 2 years apart, it was reported that stable results were obtained when the wPLI approach was applied to four 12-s epochs rather than to 4-s epochs (Hardmeier et al. 2014). In another previous study, 128-channel eyes-closed resting-state EEG data were recorded from the same participants one week apart. Connectivity indices were obtained from fixed-length 2-s epochs using the COH and wPLI approaches, and prediction stability was examined in each frequency band. The analyses found that the COH approach produced more stable results in the delta and theta bands, while the wPLI approach produced more reliable results at higher frequencies such as beta. Intra-Class Correlation (ICC) reliability was lower in the gamma band (Kuntzelman and Miskovic 2017). Based on these findings, we adopted an epoch length of 6 s, which provides a balance between capturing sufficient data points for reliable spectral and phase-based estimates. This duration is long enough to ensure stable connectivity and graph indices, but short enough to limit contamination by slow drifts, vigilance fluctuations, or artifacts that are more likely to affect longer epochs. Thus, the 6-s window represents an empirically supported compromise, consistent with previous methodological recommendations for resting-state sensor-level connectivity studies. After segmenting the data into epochs, a notch filter at 50 Hz was applied to suppress powerline interference, followed by a fourth-order Butterworth band-stop filter (60–100 Hz) to attenuate high-frequency EMG activity. And a band-pass filter (FIR filter) was applied sequentially. The band-pass filter's frequency range was adjusted to the targeted EEG sub-band. Other filter parameters were standardized across all frequency bands (pass-band ripple: 0.057501127785, stop-band attenuation: 0.0001, density factor: 20).

In the second step, two distinct approaches were employed to remove noise. In the first approach, the Artifact Subspace Reconstruction (ASR) algorithm suggested in reference (Mullen et al. 2015) was applied, while in the second, Independent Component Analysis (ICA) suggested in reference (Castellanos and Makarov 2006) was used as described in Sect. 2.3. For both approaches, the filtered and artifact-cleaned EEG segments were subjected to connectivity methods (DTF, iCOH, wPLI), which are introduced in the subsequent subsections. The DTF method was implemented using the eConnectome toolbox (Huang et al. 2021), while EEG segment modeling was performed with the AR-Fit toolbox (Schneider and Neumaier, 2001). Connectivity indices were calculated from the connectivity values computed for each 2-s segment of the 20-channel surface EEG recording using the Brain Connectivity Toolbox (Rubinov and Sporns 2010. Each method’s implementation steps are detailed and introduced in the following subsections. All computational steps were performed using MATLAB 2025Ra. First, to identify which connectivity estimation method yielded the most stable estimates across repeated epochs, we calculated Intraclass Correlation Coefficients (ICCs) for each frequency-band-specific connectivity index separately for the EO and EC states across 40 non-overlapping 6-s epochs for each participant. After determining the method that provided the most stable indices, we subsequently examined group-level statistical differences in the connectivity indices derived from that method. ICCs were computed using a MATLAB implementation of the Shrout and Fleiss framework for assessing rater and measurement reliability (Shrout & Fleiss 1979. We employed the function provided by Brownhil, in which rows represent raters or repeated measurements and columns represent the individual targets or cases being evaluated (Brownhill 2025). Following the Shrout–Fleiss schema, the ICC model was selected as ICC(3,k), appropriate when each target is measured by the same fixed set of raters and the reliability of the mean rating across raters is of interest. The resulting ICC values were interpreted according to widely adopted thresholds, defining reliability as poor for ICC < 0.40, fair for 0.40–0.60, good for 0.60–0.75, and excellent for ICC ≥ 0.75 (Jin et al. 2011; Hardmeier et al. 2014; Kuntzelman and Miskovic 2017).

Statistical differences between the two groups were analyzed using Linear Mixed-Effects models (LMEs) in MATLAB, with a significance criterion of p < 0.05. Regarding ICC values originating from the integration of three methods (DTF, iCOH, wPLI) with ASR and ICA, LMEs were determined for the specified integration, producing estimates with good or excellent reliability (ICC ≥ 0.60). In detail, separate LMEs were constructed for each band-specific connectivity index as described in Sect. 2.5. Finally, the most meaningful estimations were visualized in figures.

EEG data acquisition and ethical approval

EEG was recorded from twenty different locations (Fp1, Fp2, F3, F4, Fz, F7, F8, C3, C4, Cz, T7, T8, P3, P4, Pz, P7, P8, O1, O2, and Oz) and electrode placements based on the international 10–20 system. The EEG was amplified by Enobio32 (NeuroElectrics, Spain) with a bandpass of 0 to 125 Hz and digitized online at 500 Hz, with the right earlobe as the reference electrode. The impedance was kept below 10 kΩ during all recordings.

Participants sat in a dimly lit room during EEG recordings. They were seated in a comfortable armchair during the rsEEG recording and instructed to remain awake, with no movement, to avoid focusing on external stimuli, and to avoid cognitive tasks such as planning during recording. Once the participants understood the instructions, a 6-min eyes-open rsEEG recording was first taken. Then, a 6-min eyes-closed rsEEG recording was taken. All procedures conducted in studies involving human participants adhered to the ethical standards set forth by the institutional and/or national research committee, as well as the 1964 Helsinki Declaration and its subsequent amendments or equivalent ethical guidelines. The study protocol received approval from the Istanbul Medipol University Non-Interventional Clinical Research Ethics Committee (Approval No: E-10840098–772-02–897). Prior to data collection, researchers provided participants with a comprehensive explanation of all experimental procedures. The informed consent form, signed by all participants, also stated that the recorded EEG data would be analyzed using computer-based methods for scientific research. Informed consent was obtained from all participants, and the aforementioned ethics committee approved the study.

Participants

19 mineworkers working underground in a single mine site and 19 control participants matched with these workers in terms of age, working duration, and education level were included in this study. Since working underground is legally prohibited for women in Türkiye, the entire study sample consists of male participants. Within the scope of the study, underground mineworkers were first selected according to predetermined inclusion and exclusion criteria; then, a control group matched for demographic characteristics was selected. Previous studies were used to determine the criteria (Çelik et al. 2025; Çelik, 2025).

Inclusion criteria for mine workers: (I) Having at least 5 years of underground mining experience, (II) Having a minimum primary school degree, (III) working the day shift during the data collection period. Exclusion criteria: (I) History of major head trauma or work accident, (II) diagnosis of schizophrenia or bipolar disorder in first-degree relatives, (III) diagnosis of any psychiatric/neurological disorder, (IV) use of pharmacological agents affecting cognitive functions (psychostimulants, antidepressants, anxiolytics, etc.) within the last year, (V) excessive alcohol consumption within the last month, (VI) history of substance use within the last six months. Participants in the control group were selected from occupations that worked shifts, including night shifts, during the study period. The socio-demographic characteristics and statistical summaries of psychometric test scores of the participant groups are listed in Table 1.

Table 1 Sociodemographic characteristics and statistical summaries of psychometric test scores of the participant groupsExamining ICA and ASR

Independent Component Analysis (ICA) was performed on each preprocessed 6-s segment using the FastICA algorithm included by FastICA toolbox v2.5 (Hyvärinen et al. 1999; Gävert et al. 2005), implemented in the NIRS-KIT toolbox. For each EEG segment, the algorithm initially decomposed the data into several independent components equal to the number of recording channels. Potential artifact-related components were then automatically identified using a kurtosis threshold of 1.25 and an amplitude threshold of 4 standard deviations. These components were excluded, and the remaining independent components were back-projected to the sensor space to reconstruct the cleaned EEG signals.

In practice, we observed that none of the channels had to be rejected, so the reconstruction always retained the complete set of components equal to the number of EEG channels. This is likely due to the preceding preprocessing steps, which effectively attenuated muscle-related high-frequency activity and other large-amplitude artifacts prior to ICA. To mitigate high-amplitude artifacts (e.g., movement and muscle activity) without discarding entire data segments, we employed ASR according to the principles implemented in the clean_rawdata EEGLAB ver. 2025.0.0 plugin. ASR was developed by Kothe, the main developer of BCILAB (Kothe and Makeig 2013). Specifically, we adapted the algorithm into an independent MATLAB function for single-channel signals. ASR works by estimating a reference covariance matrix from the continuous data, performing eigendecomposition, and monitoring short sliding windows (default: 0.5 s with 50% overlap). For each window, the variance is compared against a threshold defined as the mean eigenvalue plus cutoff × standard deviation (here, cutoff = 5). If the variance exceeds this limit, the affected window is projected back into a lower-dimensional subspace (controlled by maxdims, default 0.66), thereby suppressing transient high-amplitude artifacts while preserving the underlying neural signal structure. In our analysis, EEG recordings were first segmented into 6-s epochs and preprocessed. Then, each epoch was processed independently with the ASR function, ensuring that artifact correction was applied consistently across trials. This implementation follows the same statistical rationale as the canonical EEGLAB-based ASR approach and has been shown to effectively reduce non-neural contamination in EEG (Chang et al. 2018, 2020).

Connectivity measures

Graph theory-based functional brain connectivity analysis provides a robust framework for understanding the topological properties of brain networks through various quantitative metrics so called Modularity (Q), Global Efficiency (GE), Local Efficiency (LE), Clustering Coefficient (CC), Transitivity (T) and Assortativity (R) (Bullmore and Sporns 2012; Rubinov & Sporns 2010). Q quantifies the network's tendency to divide into modules (clusters with dense internal connections and sparse inter-module connections), where high modularity signifies functional segregation and specialized information processing. GE quantifies the overall efficiency of information transfer across the network, where high GE, associated with shorter path lengths, reflects rapid and integrated information exchange between brain regions, indicative of functional integration. LE measures the efficiency of connections among a node's neighbors, quantifying local information processing capacity and thus representing functional segregation. The CC measures the degree to which a node's neighbors are interconnected; a high clustering coefficient indicates local cliquishness and segregation, reflecting the density of local functional connections. T measures the proportion of triadic structures (triangles) in the network, reflecting local connectivity density and thus segregation. R assesses the tendency of nodes to connect with others of similar degree; positive assortativity suggests network resilience, while negative assortativity may indicate vulnerability. These metrics collectively provide critical insights into the functional organization of brain networks at both local and global levels, with applications in the study of neurological disorders and cognitive performance.

In detail, functional integration refers to the brain's ability to coordinate and combine information across distributed regions to produce unified cognitive or behavioral outcomes. This process relies on the efficient communication between distant brain areas, often facilitated by long-range connections in neural networks. Mechanistically, functional integration is supported by synchronized neural activity, as measured by EEG or fMRI, and is reflected in metrics such as GE in graph-theoretic analyses. High functional integration enables the brain to perform complex tasks, such as decision-making or perception, by integrating sensory, motor, and cognitive information. For instance, during a cognitive task, the prefrontal cortex may integrate inputs from sensory areas to guide behavior, relying on coherent oscillatory activity across frequency bands. This concept is critical in understanding how the brain achieves holistic processing despite its distributed architecture (Friston 2011).

Functional segregation describes the brain's capacity to process information in specialized, localized regions, or modules, enabling distinct functional roles. This principle is rooted in the brain's modular organization, in which specific areas, such as the visual cortex and language areas, are dedicated to particular tasks. Segregation is quantified by metrics like LE, CC, or Q in graph theory, which highlight dense local connectivity within functional units. Mechanistically, segregation arises from localized neural circuits with high within-region connectivity, often driven by short-range synaptic interactions. For example, the primary visual cortex processes specific visual features (e.g., orientation) independently before integrating them with other regions. Functional segregation ensures efficient, specialized processing, which is essential for tasks requiring distinct computational roles (Bullmore & Sporns 2012; Rubinov & Sporns 2010).

Network resilience refers to the brain's ability to maintain functionality in the face of disruptions, such as lesions, noise, or pathological conditions. This property is linked to the topological structure of brain networks, particularly their redundancy and connectivity patterns. Resilient networks often exhibit high R (where nodes of similar degree connect) or robust modular structures, allowing the brain to reroute information or compensate for damage. Mechanistically, resilience is supported by neuroplasticity and the dynamic reconfiguration of neural connections, enabling the brain to adapt to challenges such as injury or disease. For example, after a stroke, resilient brain networks may reorganize to restore function by leveraging alternative pathways. This concept is crucial for understanding recovery from brain injuries and the robustness of neural systems (Bullmore & Sporns 2012).

Directed transfer function (DTF)

DTF is a multivariate connectivity estimator based on the Multivariate Autoregressive (MVAR) model, rooted in Granger causality principles. It quantifies directed information flow between channels as a function of frequency, capturing causal interactions in brain networks. The MVAR model represents EEG signals as a linear combination of past values plus noise, transforming the time domain into the frequency domain via a Z-transform to derive a transfer matrix. For the DTF computation, each EEG epoch was modeled using a MVAR framework, and the model order (p) was optimized for every epoch. Specifically, we evaluated candidate orders within a predefined range (p = 1–30) by fitting an autoregressive model of order p to the multichannel data using the arx routine from MATLAB’s System Identification Toolbox. For each candidate order, we computed the residual covariance matrix and derived the corresponding log-likelihood of the model fit. The optimal order was then selected as the one that minimized a chosen information criterion, namely the Bayesian Information Criterion (BIC). In addition, the stability of the estimated MVAR coefficients was checked using the eigenvalue criterion, and unstable solutions were penalized. The order yielding the lowest valid criterion value across the tested range was chosen as the best model order for that epoch. DTF is calculated as the normalized magnitude of this transfer matrix, specifically the element \(_(f)\) which represents the directed influence from channel \(j\) at frequency \(f\). The formula is:

$$DTF_ \left( f \right) = \frac \left( f \right) } \right|^ }} \left| \left( f \right) } \right|^ }}$$

This normalization ensures that DTF values range between 0 and 1, reflecting the relative strength of directed connectivity. DTF accounts for multivariate interactions, distinguishing direct from indirect connections, and is particularly effective for identifying frequency-specific causal relationships. DTF is inherently robust to volume conduction because it relies on phase differences between channels, which are unaffected by the instantaneous propagation of electromagnetic fields (volume conduction). Volume conduction introduces zero-phase-lag correlations, which do not contribute to DTF estimates, as DTF is non-zero only when a phase difference exists. Simulations have demonstrated that adding a constant-phase signal (e.g., a 20 Hz sinusoid) to EEG data does not alter DTF results, confirming its insensitivity to volume conduction. This robustness eliminates the need for preprocessing steps like source projection or Independent Component Analysis (ICA), which can introduce errors if misapplied (Jung et al. 2000; Kaminski et al. 2014). However, ICA is applied to noisy EEG segments.

Imaginary Part of Coherence (iCOH)

iCOH measures functional connectivity by focusing on the imaginary component of the coherence function, which quantifies the phase relationship between two signals. Coherence is defined as the normalized cross-spectral density:

$$COH_ \left( f \right) = \frac \left( f \right)}} \left( f \right)S_ \left( f \right)} }}$$

where \(_(f)\) is the cross-spectral density between channels \(i\) and\(j\), \(_(f)\) and \(_(f)\) refer the auto-spectral densities of each channel. The imaginary part of\(_\left(f\right)\), isolates interactions with non-zero phase lags, as it is maximal when the phase difference is \(\pm \pi /2\) and zero when the phase difference is 0 or \(\pi\). iCOH thus captures true neural interactions while discarding instantaneous correlations. However, iCOH may underestimate connectivity when true interactions occur at zero or \(\pi\)-phase differences, as it ignores the real part of coherence, which may contain relevant information (Sanchez et al. 2018).

iCOH is designed to be robust against volume conduction because it exclusively uses the imaginary component of coherence, which is unaffected by zero-phase-lag effects caused by volume conduction. Since volume conduction results from the instantaneous spread of electric fields, it produces correlations with no phase lag, which are filtered out by iCOH. However, its reliance on the imaginary part may cause it to miss true interactions with zero or π-phase lags, potentially reducing sensitivity in certain scenarios. Despite this, iCOH provides a reliable estimate of functional connectivity in sensor-space EEG/MEG, as demonstrated in studies comparing it to other methods (Sanchez et al. 2018).

Weighted phase lag index (wPLI)

wPLI is an extension of the Phase Lag Index (PLI), which measures the asymmetry in the distribution of phase differences between two signals. PLI is defined as:

$$PLI_ = \left| \left( t \right)} \right)} \right|$$

where \(_(t)\) is the phase difference between signals \(i\) and \(j\) at time \(t\), sign refers the signum function. wPLI improves upon PLI by weighting phase differences by the magnitude of the imaginary part of the cross-spectral density, emphasizing larger phase lags:

$$wPLI_ = \frac} \left( } \right)} \right|sign\left( } \left( } \right)} \right)} \right|}}} \left( } \right)} \right|}}$$

This weighting reduces the influence of small phase differences near zero, enhancing robustness against noise and volume conduction. wPLI values range from 0 to 1, with higher values indicating stronger phase synchronization. The debiased wPLI estimator further improves accuracy by correcting for sample-size bias (Vinck et al. 2011).

wPLI is highly effective at mitigating volume conduction effects because it focuses on phase differences and weights them by the imaginary part of the cross-spectrum, which is insensitive to zero-phase-lag correlations. By downweighting phase differences near zero, wPLI reduces false positives caused by volume conduction or common-reference artifacts. Studies have shown that wPLI provides clearer topographic connectivity maps compared to coherence, particularly in the alpha band, due to its robustness against volume conduction and noise (Ortiz et al.2012).

Estimation of graph theoretical indices

For each participant, in a given state (EO vs EC) and within a specified EEG frequency band, inter-hemispheric connectivity levels obtained from each 6-s segment of the 20-channel surface EEG recordings, using a connectivity method (DTF, iCOH, wPLI), were organized into a connectivity matrix, as illustrated in Fig. 1.

Fig. 1The alternative text for this image may have been generated using AI.

EEG electrode placement, recording channel numbers used in the construction of the connectivity matrix, and their correspondence with matrix elements

The connectivity matrices derived from the DTF, iCOH, and wPLI measures, were converted into binary adjacency matrices prior to computing network metrics. Because DTF estimates directed (asymmetric) connections, the resulting adjacency matrix is not symmetric. A proportional threshold corresponding to 60% of the maximum connection weight was applied to derive binary adjacency matrices, ensuring the preservation of the strongest functional connections while maintaining network comparability across subjects in recognizing emotional states identified by statistical and effective connectivity estimations (Kılıç & Aydın, 2022; Aydın & Onbaşı, 2024). In the present study, we applied a minimum-spanning-tree (MST) based thresholding strategy to select the most informative connections while preserving the global structure of the network. Specifically, the weighted directed DTF matrix was first inverted so that the MATLAB minspantree function could identify the tree of strongest connections. The MST ensured that all nodes remained connected using the minimal set of highest-weight edges. We then added the strongest remaining non-MST edges until a fixed target network density (e.g., 15% of all possible edges) was reached, after which the matrix was binarized for graph-theoretic computations. For the wPLI and iCOH measures, which produce undirected (symmetric) connectivity matrices, thresholding was simpler because edge directionality was not considered. After symmetrizing the matrices, we applied a proportional threshold, retaining only the top-strength connections up to a predefined density level, followed by binarization. To achieve this, all non-MST connections were sorted in descending order of strength, and the strongest edges were added until the target density was reached. The threshold for inclusion was defined by the weakest edge that met this density criterion. In the final step, the combined MST and additional strong connections were binarized, yielding a fully connected, comparable binary adjacency matrix that retained only the strongest links. This two-tiered approach ensured that, for DTF, the network remained fully connected and directionality was preserved via MST-guided thresholding, whereas for wPLI/iCOH, a standard density-controlled threshold on the symmetric matrix was sufficient. On these binary network, six different connectivity indices (Q, GE, LE, CC, T, R) were calculated using functions from the Brain Connectivity Toolbox (BCT) (Rubinov & Sorns, 2010; 2011).

For directed networks, modularity (Q) is defined as,

$$Q = \frac\mathop \sum \limits_ \left( - \frac^ .k_^ }}} \right)\delta \left( ,c_ } \right)$$

where \(_\)​ is the adjacency matrix, \(_^\) and \(_^\)​ denote the out-degree and in-degree of nodes \(i\) and \(j\), \(m\) is the total number of edges, and \(\delta \left(_,_\right)\) equals 1 when nodes belong to the same community.

Global Efficiency (GE) is defined by,

$$GE = \frac \right)}}\mathop \sum \limits_ \frac }}$$

where \(_\)​ is the shortest directed path length from node \(i\) to node \(j\). If no directed path exists, the term is treated as zero.

Local Efficiency (LE) in directed networks is defined as the average global efficiency computed on each node’s directed neighborhood subgraph in form,

$$LE = \frac\mathop \sum \limits_ GE\left( } \right)$$

where \(_\) is the subgraph composed of neighbors of node \(i\), preserving directionality.

Clustering Coefficient (CC) is obtained as the averaged network-level coefficients calculated in form,

$$CC = \frac }}^ \left( ^ - 1} \right) - 2k_^ }}$$

where \(_\) is the number of directed triangles around node \(i\), \(_^=_^+_^\), and \(_^\) is the number of reciprocal (bidirectional) connections.

For directed grpahs, Transitivity (T) is defined as,

Assortativity (R) in directed networks measures the Pearson correlation between degrees of connected node pairs. Directed assortativity can be defined based on combinations of in-degree and out-degree (e.g., out–in degree correlation). In this study, R was computed using the directed degree correlation implemented in the BCT.

Importantly, the target density was kept identical across all participants and connectivity measures to ensure direct comparability between subjects and groups. Because graph-theoretical metrics are known to be sensitive to variations in network density, enforcing a fixed proportional density prior to binarization minimized potential density-driven bias in subsequent analyses.

Group comparisons through linear mixed effect models

Recent advances in EEG analysis and network neuroscience emphasize that averaging over trials or epochs may discard valuable within-subject variability, and instead recommend using Linear Mixed Effects models (LMEs) to properly model repeated observations. For example, in a comparative framework, the researchers illustrate how LMEs often outperform classical ANOVA in EEG analysis by better accounting for intra-subject variance and reducing false positives (Frömer et al. 2018). Moreover, other works explicitly recommend replacing subject-wise averaging with LMEs when feasible, pointing to LMEs' increased sensitivity and robustness under realistic noise conditions (e.g., unequal trial counts, missing data) compared with ANOVA approaches (Riha et al. 2020). In a more recent study, Heise et al. (2022) argue that traditional averaging of epochs can bias results when trial counts differ across participants, and demonstrate that LME yields more accurate, unbiased estimates even with missing or sparse data. Thus, by adopting LME in our analysis, we align with best practices in the neurophysiology and EEG literature, modeling both fixed group differences and random subject effects while fully utilizing the epoch-level data.

In this study, to account for the repeated-measures structure of the data (40 epochs per participant across both EO and EC states) and to improve on the simple group-wise t-tests originally reported, we fitted Linear Mixed Effects models (LMEs) for each graph-theoretic index across all frequency-band combinations and eye states. Therefore, for each LME, all single epoch estimations from both groups (Controls vs Miners) were first arranged in a long format table including the variables: Value: graph metric score of a given participant at a given epoch, Group: categorical fixed factor with two levels (Control as the reference level and Miners as the second group), Subject: participant identifier to model subject-specific variability, Epoch: repeated-measurement index (1–40). The primary model specification was in form,

$$Value_ = \beta_ + \beta_ + Group_ + u_ + \varepsilon_$$

where \(_\) captures the mean difference between groups, \(_\sim N\left(0,_^\right)\) is the random intercept for \(subject i\), and \(_\sim N\left(0,^\right)\) is the residual error for \(epoch j\) within \(subject i\). LMEs were estimated using MATLAB’s fitlme function with syntax: Value ~ Group + (1 | Subject). This approach appropriately accounts for the non-independence of epoch-wise observations within each participant and relaxes the normality assumption required by ordinary t-tests.

For each state–band–index combination, the β₁ coefficient represents the estimated mean difference between the Miners group and the Control group (reference). A negative β₁ indicates lower average values in Miners compared with Control, while a positive β₁ indicates higher average values in Miners. The associated p-value tests the null hypothesis that β₁ = 0 (no group effect). This LME framework allowed us to test group differences while fully exploiting the repeated-epoch structure of the dataset and provided more reliable inference than simple independent-samples tests.