Package 'PhysioCrossModal'

Description

Establishes row correspondence before positional multiblock analyses such as jive() or multipleFactorAnalysis(). Keys may describe, for example, participant, trial and cycle together. Their meaning and granularity must be chosen by the caller; equality does not verify units, timing or independence.

Usage

alignObservationBlocks(blocks, keys, by, missing = c("error", "intersection"))

Arguments

Value

A list with blocks aligned in first-block observation order, observations containing their composite keys, row_index mapping output rows to each input block, excluded_rows listing omitted input row indices, and by and missing recording the chosen contract. Duplicate keys are always errors. The caller should retain the input keys with this mapping.

Examples

a <- matrix(1:6, nrow = 3)
b <- matrix(7:12, nrow = 3)
ids <- data.frame(participant = c("P1", "P2", "P3"))
aligned <- alignObservationBlocks(
  list(emg = a, motion = b[3:1, , drop = FALSE]),
  list(emg = ids, motion = ids[3:1, , drop = FALSE]), "participant")

Description

Takes two or more named PhysioExperiment objects and resamples them so they all share a single sampling rate, then wraps the result in a MultiPhysioExperiment.

Usage

alignSignals(
  ...,
  method = c("lowest_rate", "common_rate", "resample", "elastic"),
  target_rate = NULL
)

Arguments

Details

Three alignment strategies are available:

"lowest_rate"

(default) Resample all signals to the lowest sampling rate present.

"common_rate"

Resample all signals to the highest sampling rate present.

"resample"

Resample all signals to the rate given by target_rate (which must be supplied).

"elastic"

Resample to the lowest rate, then additionally apply an elastic (SRVF) time warping of every object onto the first, so that shared events line up in phase. The warping is computed from each object's first channel and carried to all its channels; it is stored in metadata()$elastic_warp. Objects whose duration differs from the reference after rate alignment are left un-warped. See elasticAlign().

If every input already shares the same sampling rate, no resampling is performed regardless of the chosen method (elastic warping, if requested, is still applied).

Value

A MultiPhysioExperiment containing the (possibly resampled) input objects.

See Also

alignToRate(), mergePhysio(), MultiPhysioExperiment(), couplingAnalysis()

Examples

eeg <- PhysioExperiment(
  assays = list(raw = matrix(rnorm(500 * 2), nrow = 500)),
  samplingRate = 500
)
emg <- PhysioExperiment(
  assays = list(raw = matrix(rnorm(1000 * 2), nrow = 1000)),
  samplingRate = 1000
)
mpe <- alignSignals(EEG = eeg, EMG = emg, method = "lowest_rate")

Description

Resamples all assays in a PhysioExperiment object so that the resulting data correspond to the specified target sampling rate. If the PhysioPreprocess package is available and method = "linear", its resample() function is used; otherwise an internal interpolation is applied.

Usage

alignToRate(x, target_rate, method = c("linear", "spline", "fft"))

Arguments

Value

A new PhysioExperiment with resampled data and updated samplingRate. The assay name is preserved from the original.

See Also

alignSignals(), mergePhysio(), MultiPhysioExperiment()

Examples

eeg_data <- matrix(rnorm(1000 * 4), nrow = 1000, ncol = 4)
pe <- PhysioExperiment(
  assays = list(raw = eeg_data),
  colData = S4Vectors::DataFrame(
    label = paste0("Ch", 1:4),
    type = rep("EEG", 4)
  ),
  samplingRate = 1000
)
pe500 <- alignToRate(pe, target_rate = 500)
samplingRate(pe500)  # 500

Description

Computes a bootstrap confidence interval for the coupling statistic between two signals using the moving-block bootstrap (to preserve temporal autocorrelation).

Usage

bootstrapCI(
  x,
  y = NULL,
  sr = NULL,
  method,
  n_boot = 199L,
  ci = 0.95,
  block_len = NULL,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L,
  cores = 1L,
  ...
)

Arguments

Value

A list with components:

observed

The full coupling result from the original signals.

statistic

Numeric scalar: the extracted coupling statistic.

ci_lower

Numeric scalar: lower CI bound.

ci_upper

Numeric scalar: upper CI bound.

ci_level

Numeric scalar: confidence level used.

boot_distribution

Numeric vector of bootstrap statistics.

References

Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall/CRC.

See Also

surrogateTest(), couplingAnalysis()

Examples

sr <- 500
t <- seq(0, 2, length.out = sr * 2)
x <- sin(2 * pi * 10 * t) + 0.3 * rnorm(length(t))
y <- 0.8 * sin(2 * pi * 10 * t) + 0.3 * rnorm(length(t))
result <- bootstrapCI(x, y, sr = sr, method = "coherence",
                      n_boot = 19, nperseg = 128L)
c(result$ci_lower, result$ci_upper)

Description

Finds the linear combinations of two modality blocks whose scores are maximally correlated – the directions in which the modalities co-vary. Ridge regularisation keeps it well-posed when a block has more features than observations.

Usage

cca(X, Y, ncomp = NULL, lambda_x = 0, lambda_y = 0)

Arguments

Value

a cca object: cor (canonical correlations, decreasing), xcoef, ycoef (canonical vectors), xscores, yscores (canonical variates).

References

Hotelling H (1936); Vinod HD (1976) ridge CCA.

See Also

plsBlocks(), rvCoefficient()

Examples

set.seed(1); z <- matrix(rnorm(50 * 2), 50, 2)
X <- z %*% matrix(rnorm(2 * 5), 2, 5) + matrix(rnorm(50 * 5, 0, .3), 50, 5)
Y <- z %*% matrix(rnorm(2 * 4), 2, 4) + matrix(rnorm(50 * 4, 0, .3), 50, 4)
cca(X, Y)$cor

Description

ciPLV (Bruna et al. 2018) corrects the imaginary part of the PLV for the size of its real part, giving a phase-synchrony estimate that is robust to zero-lag (volume-conduction) coupling.

Usage

ciPLV(
  x,
  y = NULL,
  sr = NULL,
  freq_band,
  order = 4L,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

A single ciPLV magnitude in ⁠[0, 1]⁠.

References

Bruna R, et al. (2018). Phase locking value revisited. J Neural Eng 15:056011.

Description

Computes the magnitude-squared coherence between two signals, which quantifies the linear frequency-domain relationship between them. The coherence is defined as:

$C_{xy}(f) = \frac{|S_{xy}(f)|^2}{S_{xx}(f) \cdot S_{yy}(f)}$

where $S_{xy}$ is the cross-spectral density and $S_{xx}$, $S_{yy}$ are the auto-spectral densities.

Usage

coherence(
  x,
  y = NULL,
  sr = NULL,
  freq_range = NULL,
  nperseg = 256L,
  noverlap = NULL,
  method = c("welch", "multitaper"),
  type = c("magnitude", "imaginary", "lagged"),
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Details

Unlike the coherence() function in PhysioAnalysis (which computes coherence across all channel pairs within a single PhysioExperiment), this function is designed for cross-modal analysis: it takes two separate signals (or PhysioExperiment / MultiPhysioExperiment objects) and computes the pairwise coherence between them.

Value

A named list with components:

coherence

Numeric vector of coherence values at each frequency: magnitude-squared coherence in $[0, 1]$ for type = "magnitude", or the (signed) imaginary / lagged coherence for the other types.

type

The coherence type computed.

frequencies

Numeric vector of corresponding frequencies in Hz.

confidence_limit

Numeric scalar giving the 95\ threshold: 1 - alpha^(1 / (L - 1)) where alpha = 0.05 and L is the number of segments (Halliday et al. 1995). Coherence values above this limit are statistically significant.

n_segments

Integer number of segments used in Welch estimation.

References

Carter, G. C. (1987). Coherence and time delay estimation. Proceedings of the IEEE, 75(2), 236–255.

Halliday, D. M., Rosenberg, J. R., Amjad, A. M., Breeze, P., Conway, B. A., & Farmer, S. F. (1995). A framework for the analysis of mixed time series/point process data – theory and application to the study of physiological tremor, single motor unit discharges and electromyograms. Progress in Biophysics and Molecular Biology, 64(2–3), 237–278.

Nolte, G., et al. (2004). Identifying true brain interaction from EEG data using the imaginary part of coherency. Clinical Neurophysiology, 115(10), 2292–2307. doi:10.1016/j.clinph.2004.04.029

Pascual-Marqui, R. D. (2007). Instantaneous and lagged measurements of linear and nonlinear dependence between groups of multivariate time series. arXiv:0711.1455.

See Also

crossSpectrum(), multitaperCoherence(), laggedCoherence(), coherenceMatrix(), couplingAnalysis()

Examples

# Two coupled 20 Hz sinusoids
sr <- 500
t <- seq(0, 10, length.out = sr * 10)
x <- sin(2 * pi * 20 * t) + 0.2 * rnorm(length(t))
y <- 0.8 * sin(2 * pi * 20 * t) + 0.2 * rnorm(length(t))
result <- coherence(x, y, sr = sr)
plot(result$frequencies, result$coherence, type = "l",
     xlab = "Frequency (Hz)", ylab = "Coherence")
abline(h = result$confidence_limit, lty = 2, col = "red")

Description

Computes magnitude-squared coherence for all pairs of channels between two sets of signals (or between all channels within a single set). Returns a matrix of peak coherence values plus the full spectra.

Usage

coherenceMatrix(
  x,
  y = NULL,
  sr = NULL,
  channels_x = NULL,
  channels_y = NULL,
  modality_x = NULL,
  modality_y = NULL,
  ...
)

Arguments

Value

A list with components:

matrix

Numeric matrix of peak coherence values, with dimensions [n_channels_x x n_channels_y].

spectra

List of lists containing the full coherence result for each pair.

frequencies

Numeric vector of frequencies from the first pair.

channel_names_x

Character vector of x channel names.

channel_names_y

Character vector of y channel names.

References

Carter, G. C. (1987). Coherence and time delay estimation. Proceedings of the IEEE, 75(2), 236–255.

See Also

coherence(), couplingMatrix(), plotCouplingMatrix(), surrogateMatrixTest()

Examples

sr <- 200; n <- sr * 2
pe1 <- PhysioExperiment(
  assays = list(raw = matrix(rnorm(n * 2), nrow = n)),
  samplingRate = sr
)
pe2 <- PhysioExperiment(
  assays = list(raw = matrix(rnorm(n * 2), nrow = n)),
  samplingRate = sr
)
result <- coherenceMatrix(pe1, pe2, nperseg = 64L)
result$matrix

Description

Computes phase-amplitude coupling over a grid of phase frequencies and amplitude frequencies, yielding a comodulogram matrix whose peak locates the dominant coupling.

Usage

comodulogram(
  x,
  y = NULL,
  sr = NULL,
  phase_freqs,
  amp_freqs,
  method = c("tort", "canolty", "ozkurt", "plv"),
  phase_bw = 2,
  amp_bw = 10,
  n_bins = 18L,
  order = 4L,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

A list with the matrix (rows = phase frequencies, columns = amplitude frequencies), phase_freqs, amp_freqs, method, and peak (the (phase, amplitude) frequency of the maximum and its value).

See Also

phaseAmplitudeCoupling(), plotComodulogram()

Examples

## Not run: 
set.seed(1); t <- (0:4999) / 500; ph <- 2 * pi * 6 * t
sig <- sin(ph) + (0.1 + (1 + cos(ph - pi)) / 2) * sin(2 * pi * 50 * t)
cm <- comodulogram(sig, sr = 500, phase_freqs = c(4, 6, 8),
                   amp_freqs = c(30, 50, 70))
cm$peak

## End(Not run)

Description

Factorises two non-negative blocks that share the observation dimension into a SHARED factor matrix W and modality-specific loadings Hx, Hy, minimising ⁠||X - W Hx||^2 + ||Y - W Hy||^2⁠. The shared W is the common latent drive recovered jointly from both modalities.

Usage

coupledNMF(
  X,
  Y,
  rank,
  max_iter = 500L,
  tol = 1e-07,
  n_restart = 5L,
  seed = NULL
)

Arguments

Value

a coupled_nmf object: shared (⁠n x rank⁠ shared factor W), Hx, Hy (⁠rank x p⁠ and ⁠rank x q⁠ loadings), vaf_x, vaf_y, vaf (combined), iterations.

References

Lee & Seung (2001) NMF; Cichocki et al. (2009) coupled/joint NMF.

See Also

plsBlocks(), jive()

Examples

set.seed(1); W <- matrix(runif(50 * 2), 50, 2)      # shared drive
X <- W %*% matrix(runif(2 * 6), 2, 6); Y <- W %*% matrix(runif(2 * 4), 2, 4)
fit <- coupledNMF(X, Y, rank = 2, n_restart = 2)
fit$vaf

Description

Dispatches to the appropriate coupling function based on the method parameter. Accepts the same flexible inputs as the underlying functions: two numeric vectors (with sr), two PhysioExperiment objects, or a MultiPhysioExperiment with named modalities.

Usage

couplingAnalysis(
  x,
  y = NULL,
  mpe = NULL,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L,
  method = c("coherence", "plv", "pli", "wpli", "granger", "crosscorrelation",
    "wavelet_coherence", "wavelet_plv", "multitaper_coherence", "dtf", "pdc",
    "transferentropy", "pac"),
  sr = NULL,
  ...,
  granger_method = c("parametric", "nonparametric")
)

Arguments

Value

The result from the dispatched coupling function. See individual function documentation for details:

• "coherence": see coherence

• "plv": see phaseLockingValue

• "pli": see phaseLagIndex

• "wpli": see weightedPLI

• "granger": see grangerCausality

• "crosscorrelation": see crossCorrelation

• "wavelet_coherence": see waveletCoherence

• "wavelet_plv": see waveletPLV

References

Carter, G. C. (1987). Coherence and time delay estimation. Proceedings of the IEEE, 75(2), 236–255.

Lachaux, J.-P., Rodriguez, E., Martinerie, J., & Varela, F. J. (1999). Measuring phase synchrony in brain signals. Human Brain Mapping, 8(4), 194–208.

Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3), 424–438.

See Also

coherence(), phaseLockingValue(), grangerCausality(), crossCorrelation(), surrogateTest()

Examples

# Numeric vectors
sr <- 500
t <- seq(0, 10, length.out = sr * 10)
x <- sin(2 * pi * 20 * t) + 0.2 * rnorm(length(t))
y <- 0.8 * sin(2 * pi * 20 * t) + 0.2 * rnorm(length(t))

# Coherence
res <- couplingAnalysis(x, y, method = "coherence", sr = sr)

# Cross-correlation
res <- couplingAnalysis(x, y, method = "crosscorrelation", sr = sr)

Description

Wraps couplingAnalysis() with a per-object memoisation cache held in the couplingResults slot of a MultiPhysioExperiment. Results are keyed by a digest of the analysis arguments (method, modalities, channels, and any extra arguments), so a repeated identical request returns the stored result without recomputing.

Usage

couplingAnalysisCached(
  mpe,
  method,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L,
  ...,
  cache = TRUE
)

Arguments

Value

A list with result (the coupling result), mpe (the object with its couplingResults cache updated) and cached (TRUE when the result was served from the cache).

See Also

couplingAnalysis(), couplingResults()

Examples

eeg <- PhysioExperiment(assays = list(raw = matrix(rnorm(1000), 500, 2)),
                        samplingRate = 250)
emg <- PhysioExperiment(assays = list(raw = matrix(rnorm(1000), 500, 2)),
                        samplingRate = 250)
mpe <- MultiPhysioExperiment(EEG = eeg, EMG = emg)
out  <- couplingAnalysisCached(mpe, "coherence", "EEG", "EMG")
out2 <- couplingAnalysisCached(out$mpe, "coherence", "EEG", "EMG")
out2$cached  # TRUE

Description

Computes any coupling method for all pairs of channels between two sets of signals, extracting a scalar statistic for each pair. This is the multi-channel generalisation of couplingAnalysis.

Usage

couplingMatrix(
  x,
  y = NULL,
  sr = NULL,
  method,
  channels_x = NULL,
  channels_y = NULL,
  modality_x = NULL,
  modality_y = NULL,
  ...
)

Arguments

Value

A list with components:

matrix

Numeric matrix of coupling values.

method

Character method used.

channel_names_x

Character vector of x channel names.

channel_names_y

Character vector of y channel names.

References

Carter, G. C. (1987). Coherence and time delay estimation. Proceedings of the IEEE, 75(2), 236–255.

See Also

couplingAnalysis(), coherenceMatrix(), plotCouplingMatrix(), surrogateMatrixTest()

Examples

sr <- 200; n <- sr * 2
pe1 <- PhysioExperiment(
  assays = list(raw = matrix(rnorm(n * 2), nrow = n)),
  samplingRate = sr
)
pe2 <- PhysioExperiment(
  assays = list(raw = matrix(rnorm(n * 2), nrow = n)),
  samplingRate = sr
)
result <- couplingMatrix(pe1, pe2, method = "coherence", nperseg = 64L)
result$matrix

Description

Returns the coupling results computed for x in this session by couplingAnalysisCached.

Usage

couplingResults(x)

## S4 method for signature 'MultiPhysioExperiment'
couplingResults(x)

couplingResults(x) <- value

## S4 replacement method for signature 'MultiPhysioExperiment'
couplingResults(x) <- value

Arguments

Details

The cache is keyed on the signal values, the clock, the assay, the channels, the method and its parameters, so a result is only returned for inputs it was actually computed from. It is held for the session rather than stored inside the container: a cache carried inside a saved object outlives any evidence that it still matches its inputs.

Value

A named list of cached results (empty when nothing is cached).

See Also

couplingAnalysisCached()

Description

Fits a maximum-likelihood common-factor model to the pooled, standardised modality blocks and flags which latent factors are SHARED across modalities (load on more than one block) versus modality-specific – a latent-variable view of multi-modal integration.

Usage

crossBlockFactor(blocks, nfactors, load_threshold = 0.4)

Arguments

Value

a crossblock_factor object: loadings (variables x factors, with a block attribute), shared (logical per factor: loads on >= 2 blocks), n_shared, and the fitted factanal.

References

Spearman/Thurstone common-factor model; McDonald (1985).

See Also

jive(), multipleFactorAnalysis()

Examples

set.seed(1); f <- matrix(rnorm(80 * 1), 80, 1)
a <- f %*% matrix(rnorm(4), 1, 4) + matrix(rnorm(80 * 4, 0, .4), 80, 4)
b <- f %*% matrix(rnorm(3), 1, 3) + matrix(rnorm(80 * 3, 0, .4), 80, 3)
crossBlockFactor(list(a = a, b = b), nfactors = 1)$shared

Description

Computes the cross-correlation between two signals at various lags, measuring their time-domain coupling. Cross-correlation quantifies how similar one signal is to a time-shifted version of another. A positive peak lag indicates that y leads x (i.e., y must be shifted forward to align with x), while a negative peak lag indicates x leads y.

Usage

crossCorrelation(
  x,
  y = NULL,
  sr = NULL,
  max_lag = NULL,
  normalize = TRUE,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Details

The function accepts numeric vectors, PhysioExperiment objects, or a MultiPhysioExperiment with named modalities, using .extract_signal_pair() internally for flexible input handling.

Value

A named list with components:

correlation

Numeric vector of cross-correlation values at each lag.

lags

Integer vector of lag values in samples.

lag_seconds

Numeric vector of lag values in seconds.

peak_lag

Integer lag (in samples) at which the absolute correlation is maximised.

peak_lag_seconds

Numeric peak lag converted to seconds.

peak_correlation

Numeric cross-correlation value at the peak lag.

References

Chatfield, C. (2004). The Analysis of Time Series: An Introduction (6th ed.). Chapman & Hall/CRC.

See Also

slidingCrossCorrelation(), coherence(), couplingAnalysis()

Examples

# Simple example: y is a delayed copy of x
sr <- 500
set.seed(1)
x <- rnorm(1000)
y <- c(rep(0, 10), x[1:990])
result <- crossCorrelation(x, y, sr = sr)
result$peak_lag       # should be 10
result$peak_lag_seconds

Description

Computes the cross-spectral density (CSD) between two signals using Welch's method. The CSD captures both the magnitude and phase relationship between two signals as a function of frequency.

Usage

crossSpectrum(
  x,
  y = NULL,
  sr = NULL,
  nperseg = 256L,
  noverlap = NULL,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

A named list with components:

csd

Complex vector of cross-spectral density values.

frequencies

Numeric vector of corresponding frequencies in Hz.

magnitude

Numeric vector of CSD magnitude (Mod(csd)).

phase

Numeric vector of CSD phase in radians (Arg(csd)).

References

Carter, G. C. (1987). Coherence and time delay estimation. Proceedings of the IEEE, 75(2), 236–255.

Halliday, D. M., Rosenberg, J. R., Amjad, A. M., Breeze, P., Conway, B. A., & Farmer, S. F. (1995). A framework for the analysis of mixed time series/point process data – theory and application to the study of physiological tremor, single motor unit discharges and electromyograms. Progress in Biophysics and Molecular Biology, 64(2–3), 237–278.

See Also

coherence(), multitaperCoherence(), couplingAnalysis()

Examples

sr <- 500
t <- seq(0, 5, length.out = sr * 5)
x <- sin(2 * pi * 10 * t) + 0.1 * rnorm(length(t))
y <- cos(2 * pi * 10 * t) + 0.1 * rnorm(length(t))
result <- crossSpectrum(x, y, sr = sr)
plot(result$frequencies, result$magnitude, type = "l",
     xlab = "Frequency (Hz)", ylab = "|CSD|")

Description

Computes the cross-wavelet transform $W_{xy} = W_x W_y^*$ of two signals via the Morlet wavelet: its modulus is the common time-frequency power and its argument is the phase lead/lag between the signals. Common power is tested against an AR(1) red-noise background (Torrence & Compo 1998; Grinsted et al. 2004) by Monte-Carlo: AR(1) surrogates matched to each signal's lag-1 autocorrelation and variance are transformed identically, and the per- frequency siglvl quantile of the surrogate cross-power sets the threshold.

Usage

crossWaveletTransform(
  x,
  y = NULL,
  sr = NULL,
  frequencies = seq(1, 40, by = 1),
  n_cycles = 7,
  significance = TRUE,
  n_surrogates = 100L,
  siglvl = 0.95,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

An object of class "cross_wavelet": a list with power (⁠time x frequency⁠ cross-wavelet modulus), phase (relative phase, radians; positive = x leads y), cross (complex XWT), frequencies, times, coi (cone of influence), and — when significance — sig_level (per frequency), sig_ratio (power / level), and significant (logical mask).

References

Torrence C, Compo GP (1998). A practical guide to wavelet analysis. Bull Amer Meteor Soc 79:61-78. Grinsted A, Moore JC, Jevrejeva S (2004). Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlin Processes Geophys 11:561-566.

Examples

set.seed(1); sr <- 100; t <- seq(0, 20, by = 1 / sr)
x <- sin(2 * pi * 10 * t) + rnorm(length(t), sd = 0.5)
y <- sin(2 * pi * 10 * t - pi / 2) + rnorm(length(t), sd = 0.5)  # y lags x
xw <- crossWaveletTransform(x, y, sr = sr, frequencies = seq(4, 16, by = 1),
                            n_surrogates = 50)
mean(xw$significant)

Description

Frequency-resolved directed connectivity via the Directed Transfer Function (Kaminski & Blinowska 1991), computed from the transfer function of a multivariate autoregressive model fitted with PhysioExperiment::mvarFit(). This cross-modal implementation is numerically identical to PhysioEEG's eegDTF on the same input. Because the transfer function inverts the whole VAR system, DTF reflects both direct and indirect pathways.

Usage

directedTransferFunction(
  x,
  y = NULL,
  sr = NULL,
  order = 5L,
  freqs = NULL,
  normalized = TRUE,
  ffDTF = FALSE,
  band = NULL,
  mvar_method = "ols",
  modality_x = NULL,
  modality_y = NULL,
  channels_x = NULL,
  channels_y = NULL
)

Arguments

Value

A list with dtf (n x n x n_freqs array), the band-averaged directed matrix (rows = targets, columns = sources), frequencies, channels, and order. For a two-channel input, dtf_xy and dtf_yx give the directional summaries.

References

Kaminski, M. J., & Blinowska, K. J. (1991). A new method of the description of the information flow in the brain structures. Biological Cybernetics, 65(3), 203-210.

See Also

partialDirectedCoherence(), transferEntropy(), grangerCausality()

Examples

set.seed(1)
X <- matrix(rnorm(1500), ncol = 3)
res <- directedTransferFunction(X, sr = 100, order = 4)
dim(res$dtf)

Description

A dependence measure that captures NON-linear as well as linear cross-modal relationships and is zero if and only if the modalities are independent (Szekely et al. 2007). In [0, 1].

Usage

distanceCorrelation(X, Y, permutations = 0L, seed = NULL)

Arguments

Value

a distance_correlation object: dcor, dcov, p_value, n.

References

Szekely GJ, Rizzo ML, Bakirov NK (2007) Ann Stat 35:2769-2794.

See Also

rvCoefficient(), representationalSimilarity()

Examples

set.seed(1); x <- rnorm(60); y <- x^2 + rnorm(60, 0, 0.3)   # nonlinear
distanceCorrelation(x, y)$dcor

Description

Aligns curve y to curve x by an optimal time warping computed in the Square-Root Velocity Function representation, separating phase (the warping) from amplitude (the residual shape difference). A pure-R dynamic-programming solver is used by default; if the fdasrvf package is available it is used for higher accuracy.

Usage

elasticAlign(
  x,
  y,
  t = NULL,
  lambda = 0,
  method = c("srvf"),
  max_step = 6L,
  smooth = FALSE,
  use_fdasrvf = NA
)

Arguments

Value

A list with gamma (the warping), aligned (y warped to x), amplitude_distance, phase_distance, the SRVFs, and the solver used.

References

Srivastava, A., Wu, W., Kurtek, S., Klassen, E., & Marron, J. S. (2011). Registration of functional data using Fisher-Rao metric. arXiv:1103.3817.

See Also

srvfMean(), warpApply()

Examples

t <- seq(0, 1, length.out = 80)
x <- exp(-((t - 0.5)^2) / (2 * 0.08^2))
y <- exp(-((t - 0.65)^2) / (2 * 0.08^2))
al <- elasticAlign(x, y)
al$amplitude_distance

Description

Concatenates several modality blocks into a single design matrix for downstream (multi-view) machine learning, scaling each block so a large or high-variance modality does not dominate the fused features.

Usage

fuseBlocks(blocks, method = c("zscore", "mfa", "frobenius"))

Arguments

Value

a fused_blocks object: matrix (⁠n x sum(p_k)⁠) and block_index (column ranges per block).

See Also

multipleFactorAnalysis(), multimodalGaitIndex()

Examples

fuseBlocks(list(a = matrix(rnorm(20), 10, 2), b = matrix(rnorm(30) * 50, 10, 3)))$matrix

Description

Computes directed coupling between two physiological signals using either autoregressive prediction or nonparametric spectral factorization. Both scalar and frequency-resolved Granger causality are supported.

Usage

grangerCausality(
  x,
  y = NULL,
  sr = NULL,
  order = 5L,
  freq_range = NULL,
  method = c("parametric", "nonparametric"),
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L,
  spectral_estimator = c("multitaper", "welch"),
  n_fft = NULL,
  segment_length = NULL,
  overlap = 0.5,
  time_bandwidth = 4,
  n_tapers = NULL,
  factor_tolerance = 1e-08,
  factor_max_iterations = 1000L,
  eigen_floor = 1e-10
)

Arguments

Details

For the time-domain (parametric) method, Granger causality from x to y is defined as:

$GC_{x \to y} = \log(\sigma^2_{restricted} / \sigma^2_{unrestricted})$

where the restricted model predicts y from its own past only, and the unrestricted model predicts y from both its own past and x's past.

For spectral Granger causality (when freq_range is specified), the AR coefficients are transformed to the frequency domain via the transfer function $H(f) = A(f)^{-1}$, and the spectral GC is computed as:

$GC_{x \to y}(f) = \log(S_{yy}^{restricted}(f) / S_{yy}^{unrestricted}(f))$

The returned values are then averaged over the specified frequency band.

The parametric frequency-band path uses separate restricted and unrestricted univariate AR models and is an approximation to the full Geweke spectral decomposition. Its established output is retained for compatibility. The nonparametric path instead estimates a complete Hermitian 2x2 cross-spectral matrix with common windows/tapers, applies Wilson's minimum-phase factorization, and evaluates Geweke directional spectra. A scalar is the trapezoidal frequency mean over the selected band, with DC and Nyquist receiving endpoint half weights.

Cross-spectral eigenvalue adjustment is limited to the reported relative eigen_floor; materially indefinite matrices fail. Wilson update convergence and spectral reconstruction are checked separately. The relative reconstruction limit is min(0.05, max(0.02, 100 * factor_tolerance)), so it remains between 2% and 5%; the observed error and limit are returned in the result. Nonparametric estimation requires at least 32 paired samples and at least two spectral realizations.

Granger causality is a stationary linear predictive-direction measure. It is not proof of physical, mechanistic, or interventional causality and does not resolve unmeasured confounding, common drivers, volume conduction, instantaneous mixing, or non-stationarity.

Value

A list with components:

gc_xy

Numeric Granger causality from x to y (x drives y). A positive value indicates that x's past helps predict y beyond y's own past.

gc_yx

Numeric Granger causality from y to x (y drives x).

net_gc

Numeric net Granger causality (gc_xy - gc_yx). Positive values indicate that x drives y more than y drives x.

order

Integer AR model order used, or NA for nonparametric estimation.

method

Character estimator used.

Nonparametric results additionally contain the frequency grid, directional spectra, estimator/factorization metadata, innovation covariance, the factorization reconstruction limit, and regularization diagnostics.

References

Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3), 424-438.

Geweke, J. (1982). Measurement of linear dependence and feedback between multiple time series. Journal of the American Statistical Association, 77(378), 304-313.

Wilson, G. T. (1972). The factorization of matricial spectral densities. SIAM Journal on Applied Mathematics, 23(4), 420-426.

Dhamala, M., Rangarajan, G., & Ding, M. (2008). Analyzing information flow in brain networks with nonparametric Granger causality. NeuroImage, 41(2), 354-362.

Dhamala, M., Rangarajan, G., & Ding, M. (2008). Estimating Granger causality from Fourier and wavelet transforms of time series data. Physical Review Letters, 100, 018701.

See Also

coherence(), phaseLockingValue()

Examples

# Create directed signals: x drives y with a 5-sample lag
set.seed(42)
n <- 5000
x <- rnorm(n)
y <- numeric(n)
for (i in 6:n) y[i] <- 0.6 * x[i - 5] + 0.4 * rnorm(1)

result <- grangerCausality(x, y, sr = 500, order = 10)
result$gc_xy   # Should be positive (x drives y)
result$net_gc  # Should be positive

Description

Decomposes several modality blocks that share an observation dimension into a low-rank JOINT structure common to all blocks, low-rank INDIVIDUAL structure specific to each, and residual noise (Lock et al. 2013). The joint scores are the integrated multi-modal representation; the per-block variance split tells you how much of each modality is shared versus idiosyncratic.

Usage

jive(
  blocks,
  rank_joint = 1L,
  rank_individual = 1L,
  center = TRUE,
  scale = TRUE,
  max_iter = 1000L,
  tol = 1e-08
)

Arguments

Value

a jive object: joint and individual (lists of ⁠n x p_k⁠ matrices), joint_scores (⁠n x rank_joint⁠, the shared component scores), variance (data frame of joint / individual / residual fraction per block), ranks, iterations.

References

Lock EF, Hoadley KA, Marron JS, Nobel AB (2013) Ann Appl Stat 7:523-542.

See Also

multipleFactorAnalysis(), cca()

Examples

set.seed(1)
shared <- matrix(rnorm(40 * 2), 40, 2)                 # common subject scores
k <- shared %*% matrix(rnorm(2 * 6), 2, 6) + matrix(rnorm(40 * 6, 0, .3), 40, 6)
e <- shared %*% matrix(rnorm(2 * 5), 2, 5) + matrix(rnorm(40 * 5, 0, .3), 40, 5)
fit <- jive(list(kin = k, emg = e), rank_joint = 2, rank_individual = 1)
fit$variance

Description

Convenience wrapper for coherence with type = "lagged", the lagged (non-instantaneous) coherence of Pascual-Marqui (2007), Im(C) / sqrt(1 - Re(C)^2) where C is the complex coherency. Like the imaginary part of coherency it is insensitive to zero-lag / volume-conduction coupling.

Usage

laggedCoherence(
  x,
  y = NULL,
  sr = NULL,
  freq_range = NULL,
  nperseg = 256L,
  noverlap = NULL,
  method = c("welch", "multitaper"),
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

A named list as returned by coherence, with type = "lagged".

References

Pascual-Marqui, R. D. (2007). Instantaneous and lagged measurements of linear and nonlinear dependence between groups of multivariate time series. arXiv:0711.1455.

See Also

coherence

Examples

sr <- 200
t <- seq(0, 10, length.out = sr * 10)
x <- sin(2 * pi * 10 * t)
y <- sin(2 * pi * 10 * t + pi / 2)          # 90-degree lag
lc <- laggedCoherence(x, y, sr = sr, freq_range = c(8, 12))
max(abs(lc$coherence))

Description

Leading Eigenvector Dynamics Analysis (Cabral et al. 2017): at each time point the instantaneous phase-coherence matrix is reduced to its leading eigenvector; these are clustered over time into a small set of recurrent connectivity states. Complements state-resolved network analysis by finding the states data-drivenly rather than from external labels.

Usage

leidaStates(
  x,
  sr,
  freq_band,
  n_states = 4L,
  order = 4L,
  nstart = 10L,
  seed = NULL
)

Arguments

Value

An object of class "leida": a list with states (per-time cluster labels), centroids (states x channels leading-eigenvector centroids), occupancy (fraction of time in each state), dwell (mean run length in samples per state), and n_states.

References

Cabral J, et al. (2017). Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest. Sci Rep 7:5135.

Examples

set.seed(1); sr <- 100; n <- 2000; t <- seq_len(n) / sr
# two regimes: channels 1-2 coherent first half, 3-4 coherent second half
base <- sin(2 * pi * 10 * t)
X <- cbind(base + rnorm(n, sd = .3), base + rnorm(n, sd = .3),
           sin(2 * pi * 10 * t + 1) + rnorm(n, sd = .3),
           sin(2 * pi * 10 * t + 1) + rnorm(n, sd = .3))
res <- leidaStates(X, sr = sr, freq_band = c(8, 12), n_states = 2, seed = 1)
res$occupancy

Description

From a leidaStates() result, computes the state-to-state transition probability matrix, the switching rate, fractional occupancy, mean dwell time, and the occupancy-weighted transition entropy (predictability of the next state).

Usage

leidaTransitions(x, sr = NULL)

Arguments

Value

An object of class "leida_transitions": a list with transition (row-stochastic ⁠state x state⁠ matrix), counts (raw transition counts), switching_rate (fraction of steps that change state), switching_rate_hz (switches per second, or NA), occupancy, dwell, and entropy (mean bits of uncertainty about the next state).

Examples

set.seed(1); sr <- 100; n <- 2000; t <- seq_len(n) / sr
base <- sin(2 * pi * 10 * t)
X <- cbind(base + rnorm(n, sd = .3), base + rnorm(n, sd = .3),
           sin(2 * pi * 10 * t + 1) + rnorm(n, sd = .3),
           sin(2 * pi * 10 * t + 1) + rnorm(n, sd = .3))
res <- leidaStates(X, sr = sr, freq_band = c(8, 12), n_states = 2, seed = 1)
leidaTransitions(res, sr = sr)

Description

Evaluates model transportability by training on all-but-one sites and testing on the held-out site (LODO). Supports continuous outcomes (family = "gaussian") and binary outcomes (family = "binomial").

Usage

lodoGeneralization(
  data,
  outcome,
  site,
  features = NULL,
  family = c("gaussian", "binomial"),
  positive_class = NULL,
  threshold = 0.5,
  min_train_rows = 20L,
  scale_features = TRUE
)

Arguments

Value

A list with:

fold_metrics

Per-site metrics.

predictions

Row-level held-out predictions.

aggregate

Mean metrics across folds.

settings

Benchmark settings and feature list.

See Also

couplingAnalysis(), surrogateTest()

Examples

set.seed(1)
n <- 120
df <- data.frame(
  site_id = rep(c("A", "B", "C"), each = 40),
  f1 = rnorm(n),
  f2 = rnorm(n)
)
df$y <- 0.8 * df$f1 - 0.4 * df$f2 + rnorm(n, sd = 0.2)

res <- lodoGeneralization(
  data = df,
  outcome = "y",
  site = "site_id",
  family = "gaussian"
)
head(res$fold_metrics)

Description

Creates a pair of sinusoidal signals with a shared oscillatory component at a specified coupling frequency and controllable noise level. Useful for demonstrating and testing spectral coherence, phase synchrony, and other cross-modal coupling measures.

Usage

make_coupled_signals(
  sr1 = 500,
  sr2 = 500,
  coupling_freq = 20,
  coupling_strength = 0.8,
  noise = 0.2,
  duration = 10
)

Arguments

Details

Both signals share the same sinusoidal component at coupling_freq, scaled by coupling_strength, with additive Gaussian noise scaled by noise.

Value

A named list with components:

x

Numeric vector – first coupled signal.

y

Numeric vector – second coupled signal.

sr_x

Numeric sampling rate of signal x.

sr_y

Numeric sampling rate of signal y.

coupling_freq

Numeric coupling frequency used.

References

Carter, G. C. (1987). Coherence and time delay estimation. Proceedings of the IEEE, 75(2), 236–255.

See Also

make_eeg_emg(), make_directed_signals(), coherence(), couplingAnalysis()

Examples

signals <- make_coupled_signals(sr1 = 500, sr2 = 500,
                                 coupling_freq = 20,
                                 coupling_strength = 0.8,
                                 noise = 0.2, duration = 10)
result <- coherence(signals$x, signals$y, sr = signals$sr_x)

Description

Creates a pair of signals where x drives y with a specified lag and coupling strength. Signal x is white noise, and y is a mixture of a lagged copy of x and independent noise. This is useful for testing directed coupling measures such as Granger causality.

Usage

make_directed_signals(n = 5000, sr = 500, lag_samples = 10, coupling = 0.7)

Arguments

Details

The generating model is:

$y(t) = \text{coupling} \cdot x(t - \text{lag\_samples}) + (1 - \text{coupling}) \cdot \epsilon(t)$

where $\epsilon(t) \sim N(0, 1)$.

Value

A named list with components:

x

Numeric vector – driving signal (white noise).

y

Numeric vector – driven signal (lagged mixture).

sr

Numeric sampling rate.

lag_samples

Integer lag used.

References

Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3), 424–438.

See Also

grangerCausality(), make_coupled_signals(), crossCorrelation(), couplingAnalysis()

Examples

signals <- make_directed_signals(n = 5000, sr = 500,
                                  lag_samples = 10, coupling = 0.7)
result <- grangerCausality(signals$x, signals$y, sr = signals$sr,
                           order = 15)
result$gc_xy   # should be positive (x drives y)
result$net_gc  # should be positive

Description

Creates a pair of PhysioExperiment objects simulating simultaneous EEG and EMG recordings with corticomuscular coherence (CMC) at a specified frequency. The first channel of each modality contains a shared oscillatory component; remaining channels contain independent noise.

Usage

make_eeg_emg(
  n_sec = 10,
  n_eeg_ch = 3,
  n_emg_ch = 2,
  sr_eeg = 500,
  sr_emg = 1000,
  cmc_freq = 20
)

Arguments

Details

This function is useful for demonstrating and testing multi-channel and multi-modal coupling analyses such as coherence matrices and MultiPhysioExperiment workflows.

Value

A named list with components:

eeg

A PhysioExperiment object with EEG data (n_eeg_ch channels at sr_eeg Hz).

emg

A PhysioExperiment object with EMG data (n_emg_ch channels at sr_emg Hz).

References

Halliday, D. M., Rosenberg, J. R., Amjad, A. M., Breeze, P., Conway, B. A., & Farmer, S. F. (1995). A framework for the analysis of mixed time series/point process data – theory and application to the study of physiological tremor, single motor unit discharges and electromyograms. Progress in Biophysics and Molecular Biology, 64(2–3), 237–278.

See Also

make_coupled_signals(), make_directed_signals(), MultiPhysioExperiment(), coherenceMatrix()

Examples

data <- make_eeg_emg(n_sec = 5, n_eeg_ch = 3, n_emg_ch = 2)
mpe <- MultiPhysioExperiment(EEG = data$eeg, EMG = data$emg)
modalities(mpe)

Description

Horizontally concatenates the channels (columns) of two PhysioExperiment objects that share the same sampling rate. Channel labels are prefixed to avoid name collisions.

Usage

mergePhysio(x, y, prefix = c("x_", "y_"))

Arguments

Details

Both objects must have the same number of time points (rows) and the same sampling rate. Only the first assay of each object is merged.

Value

A single PhysioExperiment with the columns of both inputs.

See Also

alignToRate(), alignSignals(), MultiPhysioExperiment()

Examples

pe1 <- PhysioExperiment(
  assays = list(raw = matrix(rnorm(100 * 4), nrow = 100)),
  colData = S4Vectors::DataFrame(label = paste0("A", 1:4)),
  samplingRate = 500
)
pe2 <- PhysioExperiment(
  assays = list(raw = matrix(rnorm(100 * 4), nrow = 100)),
  colData = S4Vectors::DataFrame(label = paste0("B", 1:4)),
  samplingRate = 500
)
merged <- mergePhysio(pe1, pe2)

Description

Computes Tort, Canolty, or Ozkurt phase-amplitude coupling over a declared frequency grid. Each surrogate amplitude-source signal produces one complete matrix, so cell-wise p-values share a coherent null replicate. P-values use the conservative Phipson-Smyth correction and are adjusted over the complete grid.

Usage

modulationIndex(
  x,
  phase_freqs,
  amp_freqs,
  sr = NULL,
  y = NULL,
  method = c("tort", "canolty", "ozkurt"),
  phase_bw = 2,
  amp_bw = 10,
  n_bins = 18L,
  order = 4L,
  n_surrogates = 199L,
  surrogate_type = c("timeshift", "phase", "iaaft", "aaft"),
  alpha = 0.05,
  p_adjust_method = "BH",
  seed = NULL,
  cores = 1L,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Details

PAC is an association and can be produced by waveform shape, filter leakage, edge transients, non-stationarity, or common inputs. A high or significant value does not establish cross-frequency communication or a biological mechanism. The returned threshold is the unadjusted cell-wise null quantile; the significant mask is instead defined from adjusted p-values. Benjamini-Hochberg controls the false discovery rate, not strong family-wise error. Requested and exact realized pass bands are both retained; a structured diagnostic records any adjustment at the filter's numerical boundary.

Value

A pac_comodulogram list containing the observed matrix, grid and realized bands, peak, cell-wise p-values, adjusted p-values, significance mask, null thresholds and summaries, surrogate settings, and structured diagnostic strings. Inferential matrices are NULL when n_surrogates = 0.

References

Tort, A. B. L., et al. (2010). Journal of Neurophysiology, 104, 1195-1210.

Canolty, R. T., et al. (2006). Science, 313, 1626-1628.

Ozkurt, T. E., & Schnitzler, A. (2011). Journal of Neuroscience Methods, 201, 438-443.

Phipson, B., & Smyth, G. K. (2010). Statistical Applications in Genetics and Molecular Biology, 9, Article 39.

See Also

comodulogram(), phaseAmplitudeCoupling(), plotComodulogram()

Examples

## Not run: 
set.seed(1)
t <- (0:4999) / 500
ph <- 2 * pi * 6 * t
sig <- sin(ph) + (0.1 + (1 + cos(ph - pi)) / 2) * sin(2 * pi * 50 * t)
mi <- modulationIndex(sig, c(4, 6, 8), c(40, 50, 60), sr = 500,
                      n_surrogates = 99, seed = 1)
plotComodulogram(mi)

## End(Not run)

Description

A Gait-Deviation-Index-style score (Schwartz & Rozumalski 2008) generalised across modalities: how far each subject's multi-modal gait feature vector lies from a normative reference cloud, scaled so a typical control scores ~100 and every 10 points is one standard deviation of deviation. Optionally broken down per modality to show which modality drives the deviation.

Usage

multimodalGaitIndex(
  features,
  reference,
  blocks = NULL,
  ve = 0.95,
  mean_index = 100,
  sd_index = 10
)

Arguments

Value

a multimodal_index object: index (per subject; ~100 = typical), distance, reference_mean_index, and (if blocks given) by_modality (subject x modality index matrix).

References

Schwartz MH, Rozumalski A (2008) Gait Posture 28:351-357.

See Also

jive(), fuseBlocks()

Examples

set.seed(1); ctrl <- matrix(rnorm(60 * 8), 60, 8)
subj <- rbind(matrix(rnorm(5 * 8), 5, 8),               # typical
              matrix(rnorm(5 * 8, 3), 5, 8))            # deviant
multimodalGaitIndex(subj, ctrl)$index

Description

PhysioCrossModal used to define its own MultiPhysioExperiment class, holding experiments / alignment / sampleMap / couplingResults. That class no longer exists here: the ecosystem has a single multimodal container, defined in PhysioExperiment, and this package is an analysis package that takes it as input.

Details

The constructor and every accessor keep their names and return types, so existing code runs unchanged: MultiPhysioExperiment(), experiments(), modalities(), samplingRates(), nModalities(), alignment(), alignment<-, length(), names(), [, [[. alignment() is now a view computed from the container's shared clock rather than a separately stored table, so the alignment and the signals cannot drift apart.

Containers saved under the old class are converted by migrateContainer or read with readPhysioRDS.

See Also

MultiPhysioExperiment, couplingAnalysis, couplingResults

Description

length(), names(), [ and [[ for the multimodal container are defined once, in PhysioExperiment, and are imported from there. They are not redefined here: two definitions of the same method on the same class is exactly the divergence this package's own container was unified away to prevent.

Details

One behaviour changed in the move, deliberately. The [ defined here selected a time window by computing floor(t * rate) + 1 from each modality's sampling rate alone. That ignored each stream's start offset on the shared clock, so a window taken from modalities that did not start together returned mismatched spans; it also rebuilt the alignment table from defaults, resetting every offset to zero, and discarded the sample map. Selection is now made on each sample's actual time, the clock is preserved, and each stream's offset is updated to its first retained sample. For windows whose bounds land on sample times - including all simultaneously started data with grid-aligned bounds - the selection is unchanged.

The cache accessors couplingResults() and couplingResults<- are defined in coupling-wrapper.R, where the cache itself lives.

See Also

timeWindow, MultiPhysioExperiment

Description

A global PCA over several blocks in which each block is first scaled by its leading singular value, so a large or high-variance modality cannot dominate the shared factors (Escofier & Pages 1994). Gives global factor scores on the shared observations, with each block's contribution to each factor.

Usage

multipleFactorAnalysis(
  blocks,
  ncomp = 2L,
  center = TRUE,
  scale_columns = FALSE
)

Arguments

Value

an mfa object: scores (⁠n x ncomp⁠ global factor scores), eigenvalues, explained (proportion of variance per factor), contributions (block x factor, each block's share of a factor's inertia), block_weights.

References

Escofier B, Pages J (1994) Comput Stat Data Anal 18:121-140.

See Also

jive(), rvCoefficient()

Examples

set.seed(1)
s <- matrix(rnorm(30 * 2), 30, 2)
b1 <- s %*% matrix(rnorm(2 * 5), 2, 5) + matrix(rnorm(30 * 5, 0, .3), 30, 5)
b2 <- s %*% matrix(rnorm(2 * 4), 2, 4) + matrix(rnorm(30 * 4, 0, .3), 30, 4)
multipleFactorAnalysis(list(kin = b1, grf = b2))$explained

Description

Computes magnitude-squared coherence using multitaper spectral estimation with Discrete Prolate Spheroidal Sequences (DPSS / Slepian tapers). Multitaper estimation provides lower variance than Welch's method for short signals or when frequency resolution is critical.

Usage

multitaperCoherence(
  x,
  y = NULL,
  sr = NULL,
  nw = 4,
  k = NULL,
  freq_range = NULL,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L,
  ...
)

Arguments

Value

A list with components:

coherence

Numeric vector of magnitude-squared coherence values in $[0, 1]$.

frequencies

Numeric vector of corresponding frequencies in Hz.

confidence_limit

Numeric scalar giving the 95\ threshold based on the number of tapers.

References

Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70(9), 1055–1096.

Carter, G. C. (1987). Coherence and time delay estimation. Proceedings of the IEEE, 75(2), 236–255.

See Also

coherence(), crossSpectrum(), couplingAnalysis()

Examples

sr <- 500
t <- seq(0, 10, length.out = sr * 10)
x <- sin(2 * pi * 20 * t) + 0.2 * rnorm(length(t))
y <- 0.8 * sin(2 * pi * 20 * t) + 0.2 * rnorm(length(t))
result <- multitaperCoherence(x, y, sr = sr)

Description

Computes the mutual information (Shannon, 1948) between two signals – the amount of information one signal carries about the other, a symmetric, undirected, nonlinear measure of statistical dependence. It is the information-theoretic complement of transferEntropy (directed): $I(X;Y) = \sum p(x,y)\,\log\frac{p(x,y)}{p(x)\,p(y)}$, zero iff the two signals are independent.

Usage

mutualInformation(x, y, bins = 16L)

Arguments

Details

This is the histogram estimator: each signal is discretized into bins equal-width bins and the mutual information of the resulting joint distribution is returned, in nats (natural log), matching sklearn.metrics.mutual_info_score on the same binning. The histogram estimator has a positive finite-sample bias (a shuffled surrogate has non-zero MI), so compare against a shuffled control or use more data / fewer bins for a cleaner estimate.

Value

A single numeric value, the mutual information in nats.

References

Shannon, C. E. (1948). "A mathematical theory of communication." Bell System Technical Journal, 27(3), 379–423.

See Also

transferEntropy for the directed (model-free) coupling, coherence for frequency-domain coupling.

Examples

set.seed(1)
a <- rnorm(1000); b <- a + rnorm(1000)
mutualInformation(a, b)   # > 0 (dependent)
mutualInformation(a, rnorm(1000))  # near the finite-sample floor

Description

Amplitude-envelope correlation after pairwise orthogonalisation (Hipp et al. 2012; Brookes et al. 2012): the instantaneous zero-lag (leakage) component is projected out before correlating band-limited amplitude envelopes, so shared instantaneous mixing does not create spurious coupling. Symmetrised over the two orthogonalisation directions.

Usage

orthogonalizedAEC(
  x,
  y = NULL,
  sr = NULL,
  freq_band,
  order = 4L,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

A single leakage-corrected envelope correlation in ⁠[-1, 1]⁠.

References

Hipp JF, et al. (2012). Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat Neurosci 15:884-890.

Examples

set.seed(1); sr <- 200; n <- 4000
x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)); x[is.na(x)] <- 0
y <- c(rep(0, 10), x[seq_len(n - 10)]) + rnorm(n, sd = 0.5)
orthogonalizedAEC(x, y, sr = sr, freq_band = c(8, 12))

Description

The bias-free phase-synchrony estimator of Vinck et al. (2010): $PPC = (|\sum e^{i\Delta\phi}|^2 - N) / (N(N-1))$. Unlike the PLV it is not inflated by the sample size.

Usage

pairwisePhaseConsistency(
  x,
  y = NULL,
  sr = NULL,
  freq_band,
  order = 4L,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

A single PPC value (approaches 1 for perfect locking, 0 for none).

References

Vinck M, et al. (2010). The pairwise phase consistency. NeuroImage 51:112-122.

Description

Frequency-resolved directed connectivity via Partial Directed Coherence (Baccala & Sameshima 2001), computed from the coefficient matrix of an MVAR model (PhysioExperiment::mvarFit()); numerically identical to PhysioEEG's eegPDC. Unlike DTF, PDC reflects only direct influences, so a purely indirect pathway gives PDC near zero.

Usage

partialDirectedCoherence(
  x,
  y = NULL,
  sr = NULL,
  order = 5L,
  freqs = NULL,
  generalized = TRUE,
  band = NULL,
  mvar_method = "ols",
  modality_x = NULL,
  modality_y = NULL,
  channels_x = NULL,
  channels_y = NULL
)

Arguments

Value

A list as in directedTransferFunction() with a pdc array and band-averaged matrix; for a two-channel input, pdc_xy and pdc_yx.

References

Baccala, L. A., & Sameshima, K. (2001). Partial directed coherence: a new concept in neural structure determination. Biological Cybernetics, 84(6), 463-474.

See Also

directedTransferFunction(), transferEntropy()

Examples

set.seed(1)
X <- matrix(rnorm(1500), ncol = 3)
res <- partialDirectedCoherence(X, sr = 100, order = 4)
dim(res$pdc)

Description

Quantifies how the amplitude of a fast oscillation (in amp_band) is modulated by the phase of a slow oscillation (in phase_band). The phase is taken from x and the amplitude envelope from y (or from x itself when y is NULL), each via a bandpass filter and the Hilbert transform. Supported measures are the Tort modulation index (a normalized Kullback-Leibler distance of the phase-amplitude distribution from uniform, in 0..1), the Canolty mean vector length (raw), the Ozkurt normalized mean vector length (in 0..1), and a phase-locking value.

Usage

phaseAmplitudeCoupling(
  x,
  y = NULL,
  sr = NULL,
  phase_band = c(4, 8),
  amp_band = c(30, 80),
  method = c("tort", "canolty", "ozkurt", "plv"),
  n_bins = 18L,
  order = 4L,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

A list with pac (the coupling value), method, phase_band, amp_band, and, for the Tort method, the per-bin normalized amplitude distribution and bin_centers.

References

Tort, A. B. L., et al. (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. Journal of Neurophysiology, 104(2), 1195-1210.

Canolty, R. T., et al. (2006). High gamma power is phase-locked to theta oscillations in human neocortex. Science, 313(5793), 1626-1628.

See Also

comodulogram(), plotComodulogram(), surrogateTest()

Examples

set.seed(1)
t <- (0:4999) / 500
ph <- 2 * pi * 6 * t
sig <- sin(ph) + (0.1 + (1 + cos(ph - pi)) / 2) * sin(2 * pi * 50 * t)
phaseAmplitudeCoupling(sig, sr = 500, phase_band = c(4, 8),
                       amp_band = c(35, 65))$pac

Description

Computes the Phase Lag Index between two signals. PLI measures the asymmetry of the distribution of phase differences and is insensitive to zero-lag (volume conduction) effects.

Usage

phaseLagIndex(
  x,
  y = NULL,
  sr = NULL,
  freq_band,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Details

$\text{PLI} = \left|\frac{1}{N} \sum_{t=1}^{N} \text{sign}(\sin(\phi_x(t) - \phi_y(t)))\right|$

A PLI of 0 indicates either no coupling or symmetric (zero-lag) coupling. A PLI of 1 indicates perfect non-zero-lag phase locking.

Value

A named list with components:

pli

Numeric scalar in [0, 1] – the Phase Lag Index.

phase_diff

Numeric vector of instantaneous phase differences (radians).

References

Stam, C. J., Nolte, G., & Daffertshofer, A. (2007). Phase lag index: Assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Human Brain Mapping, 28(11), 1178–1193.

See Also

phaseLockingValue(), weightedPLI(), couplingAnalysis()

Examples

sr <- 500
t <- seq(0, 2, length.out = sr * 2)
x <- sin(2 * pi * 10 * t)
y <- sin(2 * pi * 10 * t + pi / 3)
result <- phaseLagIndex(x, y, sr = sr, freq_band = c(8, 12))
result$pli

Description

Computes the Phase Locking Value between two signals, measuring the consistency of phase difference across time. A PLV of 1 indicates perfect phase locking; a PLV of 0 indicates no consistent phase relationship.

Usage

phaseLockingValue(
  x,
  y = NULL,
  sr = NULL,
  freq_band,
  method = c("hilbert", "wavelet"),
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Details

The signals are bandpass-filtered to freq_band and then Hilbert- transformed to extract instantaneous phase. PLV is computed as:

$\text{PLV} = \left|\frac{1}{N} \sum_{t=1}^{N} e^{i(\phi_x(t) - \phi_y(t))}\right|$

Value

A named list with components:

plv

Numeric scalar in [0, 1] – the Phase Locking Value.

phase_diff

Numeric vector of instantaneous phase differences (radians, [-pi, pi]).

References

Lachaux, J.-P., Rodriguez, E., Martinerie, J., & Varela, F. J. (1999). Measuring phase synchrony in brain signals. Human Brain Mapping, 8(4), 194–208.

See Also

phaseLagIndex(), weightedPLI(), waveletPLV(), couplingAnalysis()

Examples

# Two phase-locked sinusoids
sr <- 500
t <- seq(0, 2, length.out = sr * 2)
x <- sin(2 * pi * 10 * t)
y <- sin(2 * pi * 10 * t + pi / 4)
result <- phaseLockingValue(x, y, sr = sr, freq_band = c(8, 12))
result$plv

Description

The Phase Slope Index (Nolte et al. 2008) estimates the direction of information flow from the slope of the coherency phase across a frequency band, and is insensitive to instantaneous (zero-lag) mixing. Positive values indicate that x leads (drives) y.

Usage

phaseSlopeIndex(
  x,
  y = NULL,
  sr = NULL,
  freq_band,
  nperseg = 256L,
  noverlap = NULL,
  normalize = TRUE,
  modality_x = NULL,
  modality_y = NULL,
  channels_x = 1L,
  channels_y = 1L
)

Arguments

Value

A single Phase Slope Index value.

References

Nolte G, et al. (2008). Robustly estimating the flow direction of information in complex physical systems. Phys Rev Lett 100:234101. doi:10.1103/PhysRevLett.100.234101

Examples

set.seed(1); sr <- 200; n <- 4000
x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)); x[is.na(x)] <- 0
y <- c(rep(0, 8), x[seq_len(n - 8)])          # y lags x
phaseSlopeIndex(x, y, sr = sr, freq_band = c(1, 40))

Description

Displays the coherence as a function of frequency from the result of coherence. Optionally draws the 95\ a horizontal dashed line.

Usage

plotCoherenceSpectrum(
  result,
  show_threshold = TRUE,
  title = "Coherence Spectrum",
  ...
)

Arguments

Value

A ggplot object.

References

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. doi:10.1007/978-3-319-24277-4

See Also

coherence(), multitaperCoherence(), plotCouplingMatrix()

Examples

sr <- 500
t <- seq(0, 10, length.out = sr * 10)
x <- sin(2 * pi * 20 * t) + 0.2 * rnorm(length(t))
y <- 0.8 * sin(2 * pi * 20 * t) + 0.2 * rnorm(length(t))
res <- coherence(x, y, sr = sr)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plotCoherenceSpectrum(res)
}

Description

Renders a comodulogram (from comodulogram()) as a heatmap.

Usage

plotComodulogram(
  result,
  title = "Comodulogram",
  mask = NULL,
  nonsignificant_alpha = 0.2
)

Arguments

Value

A ggplot object.

See Also

comodulogram(), modulationIndex()

Examples

## Not run: 
plotComodulogram(comodulogram(sig, sr = 500,
  phase_freqs = c(4, 6, 8), amp_freqs = c(30, 50, 70)))

## End(Not run)

Description

Displays a matrix of coupling values (e.g., coherence between all channel pairs) as a colour-coded heatmap using ggplot2::geom_tile().

Usage

plotCouplingMatrix(
  mat,
  title = "Coupling Matrix",
  low_colour = "white",
  high_colour = PhysioExperiment::physioPalette(2, "sequential")[1],
  ...
)

Arguments

Value

A ggplot object.

References

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. doi:10.1007/978-3-319-24277-4

See Also

coherenceMatrix(), couplingMatrix(), surrogateMatrixTest()

Examples

mat <- matrix(runif(16), nrow = 4,
              dimnames = list(paste0("Ch", 1:4), paste0("Ch", 1:4)))
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plotCouplingMatrix(mat)
}

Description

Displays the peak cross-correlation over time from the result of slidingCrossCorrelation.

Usage

plotCouplingTimecourse(result, title = "Coupling Time Course", ...)

Arguments

Value

A ggplot object.

References

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. doi:10.1007/978-3-319-24277-4

See Also

slidingCrossCorrelation(), crossCorrelation(), plotCoherenceSpectrum()

Examples

sr <- 500
set.seed(1)
x <- rnorm(5000)
y <- c(rep(0, 10), x[1:4990])
res <- slidingCrossCorrelation(x, y, sr = sr,
                                window_sec = 1, step_sec = 0.5)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plotCouplingTimecourse(res)
}