[1910.11044] Torus Graphs for Multivariate Phase Coupling Analysis
This result agrees with the notion that hypothesis testing based on PLV, even when corrected for multiple comparisons, cannot be reliably used to control FPR for multivariate graphs because PLV measures both direct as well as indirect connectivity and thus tends to overestimate connectivity

Abstract: Angular measurements are often modeled as circular random variables, where
there are natural circular analogues of moments, including correlation. Because
a product of circles is a torus, a d-dimensional vector of circular random
variables lies on a d-dimensional torus. For such vectors we present here a
class of graphical models, which we call torus graphs, based on the full
exponential family with pairwise interactions. The topological distinction
between a torus and Euclidean space has several important consequences.
Our development was motivated by the problem of identifying phase coupling
among oscillatory signals recorded from multiple electrodes in the brain:
oscillatory phases across electrodes might tend to advance or recede together,
indicating coordination across brain areas. The data analyzed here consisted of
24 phase angles measured repeatedly across 840 experimental trials
(replications) during a memory task, where the electrodes were in 4 distinct
brain regions, all known to be active while memories are being stored or
retrieved. In realistic numerical simulations, we found that a standard
pairwise assessment, known as phase locking value, is unable to describe
multivariate phase interactions, but that torus graphs can accurately identify
conditional associations. Torus graphs generalize several more restrictive
approaches that have appeared in various scientific literatures, and produced
intuitive results in the data we analyzed. Torus graphs thus unify multivariate
analysis of circular data and present fertile territory for future research.

Figure 3: Rectangular coordinates are unable to accurately represent strong positive association between two circular random variables. (A) Scatter plot of simulated observations from a pair of dependent circular variables in rectangular coordinates, with three observations highlighted in red, black, and blue (simulated plot is similar to real data plots, but somewhat more concentrated for visual clarity; see Figure ??). The highlighted observations are shown on circles at the top of the figure (one angle as a dashed line, one as a solid line; positive dependence implies a consistent offset between the angles). While the black and blue points follow the diagonal line, the red point falls near the upper left corner due to conversion to rectangular coordinates. (B) Probability density representing the two variables, plotted both in rectangular coordinates and on the torus, with the same three points marked. On the torus, there is a single band with high probability, which wraps around and connects to itself as a Möbius strip, and all three points fall on this strip. (Introduction)Figure 4: Examples of bivariate torus graph densities and the resulting marginal distributions of phase differences upon which bivariate phase coupling measures would be based. As shown in Equation ??, the density of phase differences, p, is affected not only by coupling (through g) but also by marginal concentration (through f). As a result, bivariate phase coupling measures like PLV could give misleading results. (A) Left: bivariate torus graph density with independent angles and non-uniform marginal distributions; density is shown on the torus and flattened on [−π, π] with marginal densities on each axis. Right: analytical phase difference density (p, solid) which is a product of a direct coupling factor (g, dashed) and a marginal concentration factor (f, dotted). Here, p is concentrated solely through the marginal concentration factor f, implying bivariate phase coupling measures would indicate coupling despite the independence of X1 and X2. (B) Similar to A, but with coupling between angles and uniform marginal distributions; only in this case does p correctly reflect the coupling. (Marginal distribution of phase differences in a bivariate torus graph)Figure 5: Examples of trivariate torus graph densities and the resulting densities of phase differences for each pair of variables. As shown in Equation ??, in general, the density of phase differences, p, is affected not only by direct coupling (through g) but also by indirect connections (through h). As a result, bivariate phase coupling measures like PLV will generally reflect both direct and indirect coupling. (A) Left: ground truth graphical model, with no direct connection between X1 and X2 but an indirect connection through X3. Right: analytical phase difference densities for each pair of angles (p, solid) which are each a product of a direct coupling factor (g, dashed) and an indirect coupling factor (h, dotted). The concentration in the phase difference X1 − X2 arises solely due to indirect connections. (B) Similar to A, but with direct connection only between X1 and X2; in this case, p reflects the direct coupling. (C) Similar to A, but with direct connections between all nodes. Notice that indirect connections (h) still influence the distribution of phase differences X1 − X2 by multiplying with the direct connection term (g), which increases the concentration of p and shifts the mean (compared to g). (Marginal distribution of phase differences in a trivariate torus graph model) (Data results)Figure 6: The torus graph recovers the ground truth graph structures (top panel) from realistic simulated data sets while a bivariate phase coupling measure, phase locking value (PLV), does not (edges shown for corrected p < 0.001). Left: a 3-dimensional simulated example of cross-area phase coupling where regions X1 and X2 are not directly coupled, but are both coupled to region X3. Right: a 5-dimensional simulated example of a graph structure that could be observed for channels on a linear probe with nearest-neighbor spatial dependence. In both cases, PLV infers a fully-connected graph due to indirect connections. (Simulation results)Figure 7: In simulated data with two different underlying edge densities, the average ROC curve area under the curve (AUC) was computed across 30 simulated data sets as a function of sample size (shown along the horizontal axis). The dimension of the data is indicated by line color and different markers. Panel A demonstrates that if the true underlying graph has only 25% of all possible edges present, then even for 24 dimensional data (diamond markers), a sample size of 840 (the size of our real LFP data set) is sufficient to reach AUC above 0.9. While performance degrades when the underlying graph is more dense, panel B shows that performance is still reasonable (AUC near 0.8) for 24 dimensional data with 840 samples. (Simulation results)Figure 8: (A) Depiction of recording sites in ventrolateral prefrontal cortex (PFC) and hippocampus. (B) Preprocessing to obtain phase angles: local field potential (LFP) signals are filtered using Morlet wavelets to extract phase angles from 16 Hz oscillations at a time point of interest (two signals are shown for a single trial; repeated observations of phase angles are collected across repeated trials). (Experiment and data details)Figure 9: Torus graphs and PLV graphs from low-dimensional networks of interest in LFP data. Left: cross-region connectivity between dentate gyrus (DG), subiculum (Sub), and PFC, where the torus graph (top panel) indicated that DG and PFC are each coupled to Sub; in contrast, PLV (bottom panel) inferred a fully-connected graph. Right: within-region connectivity in CA3, where the torus graph indicated a spatial dependence structure which reflects the placement of channels along a linear probe, while PLV inferred a fully-connected graph. (Data results)Figure 10: Comparison between phase angles from three LFPs located in PFC and the theoretical torus graph distribution demonstrates that torus graphs capture the salient firstand second-order behavior present in the LFP phase angles. In contrast, the sine model fails to fit the data accurately Figure ??. (A) Along the diagonal are the marginal distributions of the phase angles. The real data are represented by blue histograms and the theoretical marginal densities from the torus graph model are overlaid as solid red traces. Two-dimensional distributions (off-diagonal) show bivariate relationships, with theoretical densities above the diagonal and real data represented using two-dimensional histograms below the diagonal. (B) Plots along the diagonal same as panel A. Below the diagonal are distributions of pairwise phase differences and above the diagonal are distributions of pairwise phase sums, represented by histograms for the real data and by solid red density plots for the theoretical torus graph model. Both the real data and theoretical distributions exhibit concentration of phase differences but not phase sums, suggesting prevalence of rotational covariance, and the sufficient statistics are very similar for the theoretical and real data. (Data results)Figure 11: Circular histograms of phase differences, in degrees, from 24-dimensional LFP data for (A) significant connections between PFC and hippocampus (specifically, CA3 and Sub) and (B) significant PFC, CA3, and Sub within-region connections, The mean phase offset with 95% confidence interval is shown as a black dot with a narrow gray interval. Observations were pooled across all significant edges within or between the regions. The within-region phase differences were tightly concentrated around zero while the PFC-hippocampus phase differences were centered below zero, indicating a possible lead-lag relationship with hippocampus leading PFC. (Summary and discussion of results)Figure 12: Bivariate torus graph densities with uniform marginal distributions shown on the torus and flattened on [−π, π] under positive, negative, or both kinds of circular covariance. (A) Coupling parameters chosen to induce only positive correlation (coupling based on phase differences). (B) Coupling parameters chosen to induce only negative correlation (coupling based on phase sums). (C) Equal amounts of positive and negative coupling result in a distribution with two isotropic modes; the superposition figure (left) is shown to provide intuition about the resulting distribution (right). (Measures of positive and negative circular dependence)Figure 13: Intuition about rotational and reflectional dependence. A) Illustration of 120◦ rotational (positive) and 120◦ reflectional (negative) dependence across three hypothetical observations. Rotational dependence implies a consistent phase offset between oscillations across trials (shaded gray angles) while reflectional dependence implies a consistent phase sum (shaded blue angles) that corresponds to a consistent phase offset between oscillations after one of the oscillations has been reflected with respect to 0◦ . B) Example of an oscillation (blue) and its reflection with respect to 0◦ (dashed red). The reflected signal is leading in time, i.e. φ = 150◦ , whereas the blue signal is lagging in time, i.e. φ = −150◦ . This demonstrates that we can think of the reflected signal as moving across time in the opposite direction. In neural data this phenomenon could arise, for example, if there was bidirectional communication. C) The lines indicate dependence between phase angles used in panel A, with rotational dependence on the left and reflectional dependence on the right and trials shown in panel A marked with stars. When phases have uniform marginal distributions across trials, but exhibit phase coupling, positive dependence is observed in the bivariate relationship on the left as a line with fixed orientation at 45◦ ; rotation produces a shift along the anti-diagonal. On the other hand, negative dependence is observed in the bivariate relationship on the right as a line with fixed orientation at 135◦ ; reflection produces a shift along the diagonal. In the current example, the respective shifts are at 120◦ ; see Figure ?? to compare with the case of 0◦ . Weaker dependence blurs the relationship line but does not change the orientation. (Measures of positive and negative circular dependence)Figure 14: In data simulated from a bivariate torus graph, the average MSE over all parameters is shown in panel (A) as a function of sample size and marginal concentration. The MSE is higher overall when marginal concentration is high. (B) ROC curves averaged over 200 simulations, half of which had no edge and half of which had an edge between the variables, as a function of marginal concentration for a fixed sample size (N=50), suggesting that structure recovery is also diminished when high marginal concentration is present. (Measures of positive and negative circular dependence)Figure 15: Further detail on the simulation results shown in Figure 2 for a sample size of 840 (matching the real LFP data). Top row: ROC curves colored by dimension for two different underlying edge densities (averaged across 30 simulations). Bottom row: averaged precision curves corresponding to the same densities as the top row. (Measures of positive and negative circular dependence)Figure 16: (A) Scatter plot of data from two channels in dentate gyrus (DG). The angles follow a pattern of positive dependence similar to the simulated data of Figure ??, which was used to demonstrate the need for circular wrapping when modeling dependent phase angles. (B) Fitted torus graph density on the plane and torus. (Measures of positive and negative circular dependence)Figure 17: Similar to Figure ??, but using the sine model as the theoretical distribution. The sine model fails to accurately fit this data set, which is evident in the bivariate dependence and sufficient statistics. (A) Along the diagonal are the marginal distributions of the three phase angles, where the real data is represented by blue histograms and the theoretical marginal densities from the sine model are overlaid as solid red traces. Two-dimensional distributions (off-diagonal) show bivariate relationships, with theoretical densities above the diagonal and real data represented using two-dimensional histograms below the diagonal. The multimodal behavior of the sine model is apparent in the two-dimensional distributions, which do not appear to match the real data. (B) Plots along the diagonal same as panel A. Below the diagonal are distributions of pairwise phase differences and above the diagonal are distributions of pairwise phase sums, represented by histograms for the real data and by solid red density plots for the theoretical torus graph model. In contrast to the torus graph model, the sine model fails to accurately capture the distributions of the sufficient statistics from these data. (Measures of positive and negative circular dependence)Figure 18: (A) Adjacency matrices for the three-dimensional LFP analysis with entries colored by p-value. PLV p-values are very small for all connections, while torus graph p-values reflect finer structure, such as PFC-Sub coupling that is apparently more salient than PFC-DG coupling. The adjacency matrix for the trivariate network is a representative combination of channels reflecting the pattern that dominates in all trivariate combinations. (B) Same as (A) but for the five-dimensional LFP analysis, where the torus graph reveals nearest-neighbor structure along the linear probe in CA3 that PLV misses. (C) Edgewise p-values for each edge in each possible trivariate graph (composed of each combination of electrodes from each of the three areas). Note that channels 6 and 7 (boxed region) of DG are on the border between DG and CA3 and may be picking up signals from CA3; omitting these channels gives stronger evidence of an overall lack of connections between DG and PFC (corresponding to the adjacency matrix in (A)). (Measures of positive and negative circular dependence)Figure 19: Adjacency matrix for PLV graph with edgewise p-values determined using Rayleigh’s test of uniformity on the circle for each pairwise phase difference. Entries are colored by p-value and, compared to the torus graph adjacency matrix Figure ??C, there is very little noticeable structure in the graph even for very small p-value thresholds (aside from a lack of edges between CA3 and DG). (Measures of positive and negative circular dependence)Figure 20: For the 24-dimensional real LFP data, the diagonal shows univariate histograms which appear to have low concentration in all cases. Below the diagonal are histograms of phase differences between pairs of angles, showing some highly concentrated distributions suggesting rotational dependence; above the diagonal are histograms of phase sums, showing very little concentration, suggesting there is not strong evidence for reflectional dependence in these data. (Measures of positive and negative circular dependence)Figure 21: Examples of simulated and real data to demonstrate the validity of the simulation process. Upper right: histograms and pairwise scatter plots, bottom left: estimated PLV matrices (color scale 0.2 to 1, with red indicating higher PLV values). (A) Simulated 5-channel data with linear probe structure. (B) Real 5-channel data from CA3. (C) Simulated 3-channel data. (D) Real 3-channel data from separate regions (DG, Sub, and PFC). (Measures of positive and negative circular dependence)Figure 22: Investigation of False Positive Rate (FPR) and False Negative Rate (FNR) for graphs of varying dimensions as sample size increases. (A) FPR and FNR for PLV (left) and torus graphs (right) using an alpha level of 0.05 for the edgewise hypothesis tests with no Bonferroni correction for the number of edges. PLV has high FPR for all sample sizes while torus graphs control the FPR; on the other hand, PLV has low FNR, but torus graphs is more conservative and for low sample sizes may be missing some edges. (B) Same as A, but with Bonferroni correction. (Measures of positive and negative circular dependence)›