---
title: "Introduction to PhysioCrossModal"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to PhysioCrossModal}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)
```

## Overview

**PhysioCrossModal** provides tools for analysing coupling, connectivity,
and synchrony *between* physiological signals of different modalities.
While single-modality analysis (e.g., EEG coherence between channels) is
handled by PhysioAnalysis, PhysioCrossModal focuses on the relationship
between distinct signal types -- for example, quantifying how cortical
EEG activity drives peripheral EMG during a motor task.

The package introduces:

- **`MultiPhysioExperiment`** -- a container that holds multiple
  `PhysioExperiment` objects recorded simultaneously at potentially
  different sampling rates.
- **Spectral coupling** -- magnitude-squared coherence and cross-spectral
  density between modalities.
- **Phase synchrony** -- Phase Locking Value (PLV), Phase Lag Index (PLI),
  and weighted PLI (wPLI).
- **Directed coupling** -- time-domain and spectral Granger causality.
- **Time-domain coupling** -- cross-correlation and sliding-window
  cross-correlation.
- A unified **`couplingAnalysis()`** wrapper that dispatches to any of
  the above methods with a single function call.

This vignette walks through a cortico-muscular coherence (CMC) analysis
as a worked example -- a common paradigm in motor neuroscience where the
goal is to quantify the functional coupling between cortical EEG and
peripheral EMG during sustained isometric contraction.

**References:** Mima & Hallett (1999), Halliday et al. (1995).

## Creating simulated EEG and EMG data

In a real experiment you would import data with `PhysioCore::readEDF()`
or similar I/O functions. Here we generate synthetic signals with a
known 20 Hz coupling (the beta-band CMC peak typically observed during
steady-force tasks).

**Note:** In real CMC research, coherence values are typically in the
range 0.1--0.4. We use a higher coupling strength here for pedagogical
clarity. Additionally, real EMG data should be rectified (full-wave
rectification) before CMC analysis.

```{r simulate-data}
library(PhysioCrossModal)

set.seed(42)

# --- Parameters ---
sr_eeg <- 500          # EEG sampling rate (Hz)
sr_emg <- 1000         # EMG sampling rate (Hz)
n_sec  <- 10           # recording duration (seconds)
cmc_freq <- 20         # coupling frequency (Hz)
cmc_strength <- 0.7    # coupling strength (0 = none, 1 = perfect)

n_eeg <- as.integer(sr_eeg * n_sec)
n_emg <- as.integer(sr_emg * n_sec)
t_eeg <- seq(0, by = 1 / sr_eeg, length.out = n_eeg)
t_emg <- seq(0, by = 1 / sr_emg, length.out = n_emg)

# Shared 20 Hz oscillation (the "cortical drive")
driver_eeg <- sin(2 * pi * cmc_freq * t_eeg)
driver_emg <- sin(2 * pi * cmc_freq * t_emg)

# Four EEG channels -- only channel 1 carries the CMC signal
n_eeg_ch <- 4
eeg_data <- matrix(rnorm(n_eeg * n_eeg_ch), nrow = n_eeg, ncol = n_eeg_ch)
eeg_data[, 1] <- eeg_data[, 1] + cmc_strength * driver_eeg

# Two EMG channels -- only channel 1 carries the CMC signal
n_emg_ch <- 2
emg_data <- matrix(rnorm(n_emg * n_emg_ch) * 0.5, nrow = n_emg, ncol = n_emg_ch)
emg_data[, 1] <- emg_data[, 1] + cmc_strength * driver_emg
```

## Building PhysioExperiment objects

Each modality is stored as a standard `PhysioExperiment` from the
PhysioCore package. The key fields are the signal matrix (time x
channels), channel metadata in `colData`, and the sampling rate.

```{r build-pe}
pe_eeg <- PhysioExperiment(
  assays   = list(raw = eeg_data),
  colData  = S4Vectors::DataFrame(
    label = paste0("EEG", seq_len(n_eeg_ch)),
    type  = rep("EEG", n_eeg_ch)
  ),
  samplingRate = sr_eeg
)

pe_emg <- PhysioExperiment(
  assays   = list(raw = emg_data),
  colData  = S4Vectors::DataFrame(
    label = paste0("EMG", seq_len(n_emg_ch)),
    type  = rep("EMG", n_emg_ch)
  ),
  samplingRate = sr_emg
)
```

## Creating a MultiPhysioExperiment

The `MultiPhysioExperiment` container bundles the two modalities
together with temporal alignment metadata. Note that the EEG and EMG
have different sampling rates (500 vs 1000 Hz) -- the container
preserves both natively.

```{r build-mpe}
mpe <- MultiPhysioExperiment(EEG = pe_eeg, EMG = pe_emg)
mpe
```

You can query the container:

```{r query-mpe}
modalities(mpe)        # "EEG" "EMG"
samplingRates(mpe)     # EEG=500 EMG=1000
nModalities(mpe)       # 2
```

## Spectral coherence (CMC)

The core analysis: magnitude-squared coherence between EEG channel 1
and EMG channel 1. Coherence quantifies the linear frequency-domain
relationship and ranges from 0 (no coupling) to 1 (perfect coupling).

Both coherence (Welch's method) and Granger causality assume weak
stationarity. For a 10-second sustained isometric contraction this is
plausible, but researchers working with dynamic tasks should verify
this assumption.

```{r coherence}
coh_result <- coherence(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  channels_x = 1L,
  channels_y = 1L,
  nperseg    = 256L
)

str(coh_result)

# Find the peak coherence near 20 Hz
peak_idx <- which.max(coh_result$coherence)
cat(sprintf("Peak coherence = %.3f at %.1f Hz\n",
            coh_result$coherence[peak_idx],
            coh_result$frequencies[peak_idx]))
```

The `confidence_limit` field provides a 95% significance threshold
based on the number of segments in the Welch estimation (Halliday et al.
1995). Coherence values above this line are unlikely to arise by chance.

You can restrict the output to a frequency range of interest:

```{r coherence-band}
coh_beta <- coherence(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  freq_range = c(13, 30)   # beta band
)
```

### Visualisation

```{r plot-coherence, fig.show="hold", eval = requireNamespace("ggplot2", quietly = TRUE)}
plotCoherenceSpectrum(coh_result, title = "EEG-EMG Coherence (CMC)")
```

## Phase Locking Value (PLV)

PLV (Lachaux et al. 1999) measures the consistency of the phase
difference between two signals within a specified frequency band. A PLV
near 1 indicates that the two signals maintain a stable phase
relationship; a PLV near 0 indicates random phase differences.

```{r plv}
plv_result <- phaseLockingValue(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  freq_band  = c(15, 25)   # beta band around 20 Hz
)

cat(sprintf("PLV = %.3f\n", plv_result$plv))
```

## Granger causality

Granger causality (Granger 1969; Geweke 1982) quantifies *directed*
coupling -- it tells you whether the past of signal X helps predict
the future of signal Y, beyond what Y's own past can predict. This is
particularly relevant for CMC because cortico-muscular drive is thought
to be predominantly (though not exclusively) descending (EEG drives EMG).

```{r granger}
gc_result <- grangerCausality(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  order      = 10L        # AR model order
)

cat(sprintf("GC (EEG -> EMG) = %.4f\n", gc_result$gc_xy))
cat(sprintf("GC (EMG -> EEG) = %.4f\n", gc_result$gc_yx))
cat(sprintf("Net GC          = %.4f\n", gc_result$net_gc))
# Positive net_gc indicates EEG drives EMG more than the reverse
```

Note: for Granger causality to work correctly, both signals must be
resampled to the same rate. When you pass a `MultiPhysioExperiment`,
the internal signal extraction handles this automatically.

## The couplingAnalysis() wrapper

If you want a uniform interface to all coupling methods, use
`couplingAnalysis()`. It dispatches to the appropriate function based
on the `method` argument:

```{r wrapper}
# Coherence via wrapper
res_coh <- couplingAnalysis(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  method     = "coherence"
)

# Cross-correlation via wrapper
res_cc <- couplingAnalysis(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  method     = "crosscorrelation"
)
cat(sprintf("Peak cross-correlation = %.3f at lag = %d samples (%.1f ms)\n",
            res_cc$peak_correlation,
            res_cc$peak_lag,
            res_cc$peak_lag_seconds * 1000))
```

The wrapper accepts the same `x`/`y` and `mpe` arguments as the
individual functions, plus `method`-specific parameters via `...`.

## Working with raw numeric vectors

All coupling functions also accept plain numeric vectors with an
explicit `sr` argument. This is useful for quick one-off analyses or
when your data is not yet stored in a PhysioExperiment:

```{r vectors}
set.seed(123)
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))

# Direct vector interface
coh <- coherence(x, y, sr = sr)
gc  <- grangerCausality(x, y, sr = sr, order = 10)
cc  <- crossCorrelation(x, y, sr = sr)
```

## Signal alignment utilities

When modalities have different sampling rates, you may need to resample
before certain analyses. PhysioCrossModal provides helper functions:

```{r alignment}
# Resample a single PE to a target rate
pe_emg_500 <- alignToRate(pe_emg, target_rate = 500)

# Align multiple PEs to the lowest rate and bundle into an MPE
mpe_aligned <- alignSignals(EEG = pe_eeg, EMG = pe_emg,
                            method = "lowest_rate")

# Merge two same-rate PEs into a single PE (column-bind)
merged_pe <- mergePhysio(pe_eeg, alignToRate(pe_emg, target_rate = 500),
                         prefix = c("EEG_", "EMG_"))
```

## Statistical significance testing

PhysioCrossModal v0.2.0 provides two complementary approaches for
assessing the statistical significance of coupling estimates:
surrogate-based null-hypothesis testing and block-bootstrap confidence
intervals.

### Surrogate testing

`surrogateTest()` generates a null distribution by randomising the
phase spectrum (or circularly shifting) one signal, thereby destroying
the cross-signal relationship while preserving autocorrelation:

```{r surrogate}
surr <- surrogateTest(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  method     = "coherence",
  n_surrogates = 19,
  nperseg = 256L
)
cat(sprintf("Observed = %.3f, p = %.4f, 95%% threshold = %.3f\n",
            surr$statistic, surr$p_value, surr$threshold_95))
```

### Bootstrap confidence intervals

`bootstrapCI()` uses the moving-block bootstrap (preserving temporal
autocorrelation) to construct confidence intervals:

```{r bootstrap}
boot <- bootstrapCI(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  method     = "coherence",
  n_boot     = 19,
  ci         = 0.95,
  nperseg    = 256L
)
cat(sprintf("95%% CI: [%.3f, %.3f]\n", boot$ci_lower, boot$ci_upper))
```

## Wavelet coherence

While Welch-based coherence gives a single frequency-domain estimate
averaged over the entire recording, wavelet coherence
(`waveletCoherence()`) reveals *time-frequency* coupling dynamics.
This is particularly useful when coupling strength varies over time
(e.g. during movement onset or task transitions).

```{r wavelet}
wcoh <- waveletCoherence(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  frequencies = seq(10, 35, by = 1),
  n_cycles = 7,
  smoothing_cycles = 3
)
# wcoh$coherence is a [time x frequency] matrix
# wcoh$phase contains the phase differences
# wcoh$coi contains the Cone of Influence frequencies

# Wavelet PLV is also available:
wplv <- waveletPLV(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  frequencies = seq(10, 35, by = 1)
)
```

### Wavelet coherence visualization

```{r plot-wavelet, fig.show="hold", eval = requireNamespace("ggplot2", quietly = TRUE)}
plotWaveletCoherence(wcoh, title = "EEG-EMG Wavelet Coherence")
```

## Multi-channel coupling matrices

When you want to compute coupling across all channel pairs (e.g., to
identify which EEG electrode shows the strongest coherence with which
EMG channel), use `coherenceMatrix()` or the generic `couplingMatrix()`:

```{r matrix}
mat <- coherenceMatrix(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  nperseg    = 256L
)
# mat$matrix is a [n_eeg_channels x n_emg_channels] matrix

# Visualise as a heatmap
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plotCouplingMatrix(mat$matrix, title = "EEG-EMG Peak Coherence")
}
```

The generic `couplingMatrix()` works with any coupling method:

```{r generic-matrix}
cc_mat <- couplingMatrix(
  mpe,
  modality_x = "EEG",
  modality_y = "EMG",
  method = "crosscorrelation"
)
```

## Summary

| Method | Function | What it measures |
|--------|----------|------------------|
| Spectral coherence | `coherence()` | Linear frequency-domain coupling [0, 1] |
| Cross-spectral density | `crossSpectrum()` | Magnitude and phase of frequency coupling |
| Phase Locking Value | `phaseLockingValue()` | Phase consistency in a frequency band [0, 1] |
| Phase Lag Index | `phaseLagIndex()` | Phase coupling robust to volume conduction [0, 1] |
| Weighted PLI | `weightedPLI()` | Noise-weighted phase coupling [0, 1] |
| Granger causality | `grangerCausality()` | Directed (causal) coupling |
| Cross-correlation | `crossCorrelation()` | Time-domain coupling with lag estimation |
| Sliding cross-correlation | `slidingCrossCorrelation()` | Time-varying coupling dynamics |
| **Wavelet coherence** | `waveletCoherence()` | Time-frequency coherence via Morlet wavelets |
| **Wavelet PLV** | `waveletPLV()` | Time-frequency phase locking |
| **Surrogate test** | `surrogateTest()` | Null-hypothesis significance testing |
| **Bootstrap CI** | `bootstrapCI()` | Confidence intervals via block bootstrap |
| **Multitaper coherence** | `multitaperCoherence()` | Coherence via DPSS tapers (low variance) |
| **Coherence matrix** | `coherenceMatrix()` | Multi-channel coherence across all pairs |
| **Coupling matrix** | `couplingMatrix()` | Multi-channel coupling (any method) |
| **Matrix significance** | `surrogateMatrixTest()` | Surrogate test with FDR correction |
| **Wavelet plot** | `plotWaveletCoherence()` | Time-frequency heatmap with COI |
| Unified wrapper | `couplingAnalysis()` | Dispatches to any of the above |

For further details on each function, see the package reference manual
(`?coherence`, `?phaseLockingValue`, etc.).
