Metric-T: checking whether a connectivity result depends on the metric

When to use Metric-T

EEG functional connectivity studies often report that one group shows “higher connectivity” than another. The direction of that difference can depend on which connectivity metric was used. The weighted phase lag index (wPLI) responds only to lagged coupling, whereas magnitude-squared coherence (MSC) also responds to zero-lag coupling such as volume conduction. Both metrics are computed from the same cross-spectrum, yet they can point in opposite directions.

Use Metric-T when you want to know:

  • whether the direction of a group difference in connectivity is the same for wPLI and MSC, or reverses between them (“directional fragility”);
  • whether a disagreement between the two metrics is larger than would be expected by chance under random group labels;
  • with small samples, exactly how unusual the observed grouping is among all possible groupings, or whether a continuous covariate (e.g. years since onset) is associated with a metric-dependent direction.

Metric-T is a diagnostic that sits next to a connectivity analysis. It does not replace a test of the connectivity difference itself.

The statistic

For each metric, direction consistency (DC) is the percentage of channel-pair features for which the mean of group 1 exceeds the mean of group 2. DC is 50% when there is no consistent direction. Metric-T is

\[T = \mathrm{DC}(\mathrm{wPLI}) - \mathrm{DC}(\mathrm{MSC}),\]

in percentage points (pp), ranging from -100 to +100. A reversal is flagged when the two DC values fall on opposite sides of 50%. Significance is assessed by aligned label permutation: the same relabelling is applied to every band and comparison, which permits Westfall-Young max-T adjustment across them.

Reference: Kimura A (2026). Metric-T: A permutation-based diagnostic for directional fragility in EEG functional connectivity analysis. Neuroscience Informatics 6:100286. doi:10.1016/j.neuri.2026.100286

library(metricT)

Input data

All analysis functions take a feature list: a named list with one element per frequency band, each holding two matrices of the same size, wpli and msc, with subjects in rows and channel pairs in columns. A group vector gives one label per subject, in the same order as the rows. Missing values are not allowed.

There are four ways to obtain it.

  1. From preprocessed EEG, one CSV per subject (samples in rows, channels in columns; a time column is ignored). Filtering, artefact removal and re-referencing are done beforehand in your EEG software.

    feats <- mt_features_from_eeg_csv(files, sfreq = 128)

  2. From epoched arrays (epochs x channels x samples per subject), with mt_epoch(), the optional mt_clean_epochs(), and mt_features_from_epochs().

  3. From a feature table in CSV with columns subject, group, and one column per feature named <band>__<metric>__<pair>, for example alpha__wpli__C3-C4 and alpha__msc__C3-C4. Read it with mt_read_features_csv(); write one with mt_write_features_csv().

  4. From features.npz produced by the published Python pipeline, with mt_read_npz() (no Python required).

The default bands are:

mt_bands()
#> $theta
#> [1] 4 8
#> 
#> $alpha
#> [1]  8 13
#> 
#> $beta
#> [1] 13 30
#> 
#> $broadband
#> [1]  1 40

In this vignette, simulated data stand in for real recordings:

sim <- mt_simulate(n = c(A = 20, B = 20), seed = 1)
names(sim$features)
#> [1] "theta"     "alpha"     "beta"      "broadband"
dim(sim$features$alpha$wpli)   # subjects x channel pairs
#> [1] 40 28
table(sim$group)
#> 
#>  A  B 
#> 20 20

The CSV route round-trips the same object:

f <- tempfile(fileext = ".csv")
mt_write_features_csv(sim$features, sim$group, f)
dat <- mt_read_features_csv(f)
names(read.csv(f, check.names = FALSE))[1:4]
#> [1] "subject"              "group"                "theta__wpli__ch1-ch2"
#> [4] "theta__wpli__ch1-ch3"
unlink(f)

Two-group analysis

res <- metric_t(dat$features, dat$group, n_perm = 2000, seed = 42)
res
#> Metric-T: T = DC(wPLI) - DC(MSC), percentage points
#> Groups: A (n=20), B (n=20)
#> Permutations: 2000 (aligned label permutation, two-sided), seed 42
#> Westfall-Young family: per band
#> 
#>  Comparison      Band DC(wPLI)% DC(MSC)%  T(pp) Reversal     p p_min p_Holm
#>      A vs B     theta      50.0     39.3  +10.7          0.363 0.000  0.727
#>      A vs B     alpha       0.0    100.0 -100.0      YES 0.001 0.001  0.006
#>      A vs B      beta      64.3     57.1   +7.1          0.594 0.000  0.727
#>      A vs B broadband       0.0    100.0 -100.0      YES 0.001 0.001  0.006
#>   q_BH  p_WY
#>  0.484 0.363
#>  0.003 0.001
#>  0.594 0.594
#>  0.003 0.001
#> 
#> Directional reversals (DC values on opposite sides of 50%): 2 of 4 conditions

In this simulation the alpha and broadband results reverse: DC(wPLI) = 0%, so group B is higher on every wPLI feature, while DC(MSC) = 100%, so group A is higher on every MSC feature. Thus a study using only one of the two metrics would report the opposite direction to a study using the other. Theta and beta show no significant metric dependence.

With more than two groups, give comparisons (a list of pairs) and optionally control; the Westfall-Young family is then each band’s comparisons against the control by default (wy_family = "control"), or each band ("band"), or all conditions together ("all").

The results table is available as a data frame:

head(as.data.frame(res)[, c("band", "DC_wPLI", "DC_MSC", "T", "reversal", "p_raw", "p_min", "p_WY")])
#>                band  DC_wPLI    DC_MSC           T reversal      p_raw
#> theta         theta 50.00000  39.28571   10.714286    FALSE 0.36331834
#> alpha         alpha  0.00000 100.00000 -100.000000     TRUE 0.00149925
#> beta           beta 64.28571  57.14286    7.142857    FALSE 0.59420290
#> broadband broadband  0.00000 100.00000 -100.000000     TRUE 0.00149925
#>                  p_min       p_WY
#> theta     0.0004997501 0.36331834
#> alpha     0.0014992504 0.00149925
#> beta      0.0004997501 0.59420290
#> broadband 0.0014992504 0.00149925
plot(res)

Reading the output

  • DC(wPLI), DC(MSC): the percentage of features in which group 1 is higher. Values near 50% mean no consistent direction for that metric.
  • T: the gap between the two DC values in pp. A large |T| means the two metrics disagree about the direction.
  • Reversal: the two metrics point in opposite directions. A reversal with a small p-value is evidence that the reported direction depends on the metric choice.
  • p, p_Holm, q_BH, p_WY: unadjusted and adjusted two-sided permutation p-values. p_WY (Westfall-Young max-T) uses the dependence between conditions and is usually the preferred adjustment.
  • p_min: the smallest p-value attainable in that condition. Because DC keeps only the sign of each mean difference, strongly correlated channel pairs can make even |T| = 100 pp non-significant. If p_min exceeds your alpha, a large p-value is not evidence that the direction is stable; the design simply cannot detect fragility there.

Small samples: exact enumeration

With few subjects, all label assignments can be enumerated, so there is no Monte Carlo error and the attainable minimum p is known exactly.

small <- mt_simulate(n = c(A = 6, B = 7), n_ch = 5, n_epochs = 12, seed = 7)
ex <- metric_t(small$features, small$group, exact = TRUE, wy_family = "all")
ex
#> Metric-T: T = DC(wPLI) - DC(MSC), percentage points
#> Groups: A (n=6), B (n=7)
#> Label assignments: all 1716 enumerated (exact, two-sided)
#> Westfall-Young family: all conditions
#> 
#>  Comparison      Band DC(wPLI)% DC(MSC)% T(pp) Reversal     p p_min p_Holm
#>      A vs B     theta      50.0     60.0 -10.0          0.810 0.001  1.000
#>      A vs B     alpha       0.0     40.0 -40.0          0.512 0.001  1.000
#>      A vs B      beta      40.0     50.0 -10.0          0.768 0.001  1.000
#>      A vs B broadband       0.0     40.0 -40.0          0.520 0.002  1.000
#>   q_BH  p_WY
#>  0.810 0.953
#>  0.810 0.677
#>  0.810 0.953
#>  0.810 0.677
#> 
#> Directional reversals (DC values on opposite sides of 50%): 0 of 4 conditions

The complete labelling space

mt_label_space() computes Metric-T for all \(2^n - 2\) two-group labellings of the subjects (by default up to n = 16). mt_locate() then places any particular split on that distribution: the observed groups, an alternative cut-point, an exclusion, or a nuisance split such as recording session.

sp <- mt_label_space(small$features)
sp
#> Complete labelling space: 13 subjects, 8190 two-group labellings
#> 
#>       band n_features  SD_T m_eff_wPLI m_eff_MSC cor_DC reversal_share
#>      theta         10 20.54      15.96      8.34   0.08           0.25
#>      alpha         10 41.03       1.62      4.83   0.21           0.37
#>       beta         10 17.04      15.30     10.34   0.29           0.18
#>  broadband         10 42.23       1.99      4.82  -0.01           0.48
#> 
#> m_eff: effective number of independent features (2500 / Var(DC)).
#> reversal_share: labellings with DC(wPLI) and DC(MSC) on opposite sides of 50%.
mt_locate(sp, small$group == "A")       # the observed grouping
#>                band   T DC_wPLI DC_MSC share_all   p_exact        p_min
#> theta         theta -10      50     60 0.8021978 0.8100233 0.0005827506
#> alpha         alpha -40       0     40 0.4979243 0.5122378 0.0005827506
#> beta           beta -10      40     50 0.7626374 0.7680653 0.0005827506
#> broadband broadband -40       0     40 0.5225885 0.5198135 0.0017482517
#>           n_same_size
#> theta            1716
#> alpha            1716
#> beta             1716
#> broadband        1716
session <- rep(c("s1", "s2", "s3"), c(4, 4, 5))
mt_locate(sp, session == "s1")          # a recording session treated as a group
#>                band   T DC_wPLI DC_MSC share_all   p_exact       p_min
#> theta         theta   0      50     50 1.0000000 1.0000000 0.001398601
#> alpha         alpha -30       0     30 0.6620269 0.6447552 0.001398601
#> beta           beta   0      60     60 1.0000000 1.0000000 0.001398601
#> broadband broadband -40       0     40 0.5225885 0.5482517 0.001398601
#>           n_same_size
#> theta             715
#> alpha             715
#> beta              715
#> broadband         715

The summary reports the effective number of independent features (2500 / Var(DC)) and the share of arbitrary labellings in which the two metrics already disagree in direction. If a recording session produces a T as large as the groups of interest, batch effects should be considered before interpreting a reversal.

Continuous covariates

metric_t_trend() defines DC as the percentage of features positively associated with a numeric covariate, and permutes the covariate only within strata (for example sessions), so that differences between strata are not treated as chance.

set.seed(11)
years <- round(runif(13, 0, 40))
tr <- metric_t_trend(small$features, years, strata = session)
tr
#> Metric-T for a continuous covariate (rank association)
#> Strata: s1 (n=4), s2 (n=4), s3 (n=5)
#> Arrangements: all 69120 enumerated within strata (exact, two-sided)
#> 
#>  covariate      band n_features DC_wPLI DC_MSC  T     p p_min p_family
#>          x     theta         10      40     20 20 0.468 0.000    0.949
#>          x     alpha         10      90     60 30 0.562 0.001    0.781
#>          x      beta         10      90     70 20 0.347 0.000    0.949
#>          x broadband         10      60     60  0 1.000 0.001    1.000
#> 
#> p_family: adjusted over all rows by the maximum |T| of each arrangement.

method = "rank" (default) is invariant to monotone transformations of the covariate and the features; method = "linear" uses within-stratum centred raw values.

Graphical interface

run_metricT_app()   # requires the 'shiny' package

The app accepts the built-in simulation, a feature file (CSV or features.npz), or one EEG CSV per subject with a participants table.

Common mistakes

  • Passing raw, unpreprocessed EEG. metricT computes wPLI and MSC but does not filter, remove artefacts, or re-reference.
  • Rows of the feature matrices in a different order from group.
  • Interpreting a non-significant result as stability when p_min is above alpha.
  • Using too few permutations for the adjusted p-values; use n_perm = 10000 or more for reporting, or exact = TRUE for small samples.

Citation

citation("metricT")