Crosssignal analyses
Command  Description 

COH 
Pairwise channel coherence 
CORREL 
Pairwise channel correlation 
CC 
Coupling (dPAC) and connectivity (wPLI) 
PSI 
Phase slope index (PSI) connectivity metric 
MI 
Mutual information 
TSYNC 
Crosscorrelation and phase delay 
COH
Calculates pairwise spectral coherence across channels
The COH
function calculates the magnitudesquared
coherence
for pairs of signals with similar sampling rates.
Parameters
Parameter  Example  Description 

sig 
sig=C3,C4 
An optional parameter, to specify which channels to calculate pairwise coherence between 
sig1 
sig1=C3,C3 
In stead of sig , specify all sig1 by sig2 channel pairs 
sig2 
sig2=F3,F4 
As above 
sr 
sr=125 
Optionally, set sample rates to sr if different for a particular channel (i.e. to ensure that all signals have similar sampling rates 
spectrum 
spectrum 
Output full crossspectra coherence rather than just in frequency bands (delta, theta, etc) 
epoch 
epoch 
A flag to output coherence statistics for per epoch as well as the whole signal 
epochspectrum 
epochspectrum 
A flag to output full crossspectra coherence statistics for each epoch/channel pair 
Warning
If requesting the full spectrum
and epoch
level analyses for a large number of channels, the output database may be very large... In this case, you should use
texttable mode output (with t
) for better performance when working with the output (as destrat can be slow othewise, for datasets with
very large numbers of strata).
Using sig1
and sig2
instead of sig
can be more efficient, if you
are only interested in particular combinations of pairs. For example,
to see the pairwsie coherences between one channel and four others:
sig1=C3 sig2=F3,F4,O1,O2
sig=C3,F3,F4,O1,O2
would
evaluate 10 pairs, even if some pairwise comparisons are not of
interest).
Output
Coherence for power bands (B
x CH1
x CH2
)
Variable  Description 

COH 
Magnitudesquared coherence 
ICOH 
Imaginary coherence 
LCOH 
Lagged coherence 
Full crossspectra coherence (option: spectrum
, strata: F
x CH1
x CH2
)
Variable  Description 

COH 
Magnitudesquared coherence 
ICOH 
Imaginary coherence 
LCOH 
Lagged coherence 
Coherence for power bands, per epoch (option: epoch
, strata: E
x B
x CH1
x CH2
)
Variable  Description 

COH 
Magnitudesquared coherence 
ICOH 
Imaginary coherence 
LCOH 
Lagged coherence 
Full crossspectra, per epoch (option: epoch
and spectrum
, strata: E
x F
x CH1
x CH2
)
Variable  Description 

COH 
Magnitudesquared coherence 
ICOH 
Imaginary coherence 
LCOH 
Lagged coherence 
nb. references for methods to be added
Example
To estimate the crossspectra (up to 20 Hz) between the two EEG channels during NREM2
sleep for the tutorial individual nsrr02
:
luna s.lst nsrr02 o out.db s "MASK ifnot=NREM2 & RE \
& COH spectrum max=20 sig=EEG,EEG(sec)"
Loading this into lunaR,
k < ldb( "out.db" )
d < lx(k,"COH","CH1","CH2","F")
plot( d$F, d$COH , type="l" ,
xlab="Frequency (Hz)" , ylab="Coherence" ,
ylim=c(0,1) , lwd=2, col="red")
Following other reports (for example), we see peaks in EEG coherence especially for the low delta and sigma frequencies during NREM sleep. Naturally, the precise interpretation of coherence analysis (like any sleep EEG analysis) depends on the exact choice of montages, referencing and other technical and subjectlevel details, not to mention the potential for artifact (e.g. as illustrated here).
CORREL
Calculates pairwise Pearson correlation between signals
CORREL
estimates pairwise correlation coefficients between signals,
either for the whole signal, or epochbyepoch. When epochlevel
statistics are requested, Luna also reports the mean and median of all
perepoch statistics for a given channel pair.
Parameters
Parameter  Example  Description 

sig 
sig=C3,C4,F3,F4 
Optionally specify which signals to correlate (otherwise, do all) 
sig1 
sig=C3,C4 
Alternatively, only evaluate all sig1 by sig2 pairs 
sig2 
sig=F3,F4 
As above 
sr 
sr=128 
Resample channels to this sample rate, if they are not already at that rate 
epoch 
epoch 
Display perepoch correlations, and estimate mean and median correlation across epochs 
Output
Wholesignal correlations for pairs of channels (strata: CH1
x CH2
)
Variable  Description 

R 
Pearson product moment correlation 
R_MEAN 
(If epoch is specified) the mean of epochlevel correlations 
R_MEDIAN 
(If epoch is specified) the median of epochlevel correlations 
Wholesignal correlations for pairs of channels (option: epoch
, strata: E
x CH1
x CH2
)
Variable  Description 

R 
Perepoch Pearson product moment correlation 
Example
Here we consider the epochbyepoch correlation between the two EOG
channels for tutorial subject nsrr02
. Using
lunaC to estimate and report correlations at
the epoch level:
luna s.lst nsrr02 o out.db s "MASK if=wake & RE \
& STAGE \
& CORREL epoch verbose sig=EOG(L),EOG(R)"
STAGE
command, which we can use later when plotting results:
Using lunaR to examine the output (in R):
library(luna)
We'll attach the output database just created:
k < ldb( "out.db" )
Examining the contents with lx()
:
lx(k)
CORREL : CH1_CH2 CH1_CH2_E
MASK : EPOCH_MASK
RE : BL
STAGE : E
First, we'll extract a simple vector of sleep stage (the STAGE
variable from the STAGE
command):
ss < lx( k , "STAGE" , "E" )$STAGE
Second, we'll extract the epochbyepoch correlations:
d < lx( k , "CORREL" , "CH1" , "CH2" , "E" )
Viewing dataframe d
, note that the epochs (E
) do not start at 1,
because we previously restricted this analysis to sleep epochs:
head(d)
ID CH1 CH2 E R
1 nsrr02 EOG(L) EOG(R) 92 0.5052592
2 nsrr02 EOG(L) EOG(R) 93 0.7998685
3 nsrr02 EOG(L) EOG(R) 98 0.4698097
4 nsrr02 EOG(L) EOG(R) 99 0.3147986
5 nsrr02 EOG(L) EOG(R) 100 0.3153311
6 nsrr02 EOG(L) EOG(R) 101 0.6581321
Finally, we can plot these epochs:
plot( d$E , d$R , col = lstgcols( ss ) , pch=20 , ylab="LR EOG correlation", xlab="Epoch", ylim=c(1,1) )
abline( h=0 , lty=2 )
For this particular individual, there seems to be a clear and interesting pattern, in which we see that REM epochs (in red) are more likely to show negative correlations between the left and right EOG, reflecting the deviations due to eyemovements during REM.
CC
Calculates coupling and connectivity metrics
The CC
command estimates crossfrequency coupling (currently
dPAC) and, for
interchannel connectivity, the weighted phase lag index
(wPLI). It implements a
timeshifting randomization to generate the null distributions of
these metrics.
Note
Although CC
currently only implements these two
metrics, in future releases other coupling and connectivity metrics
will be added within this framework (e.g. the coherence metrics
currently performed using COH
).
For one of more signals, phaseamplitude coupling is estimated by the
dPAC method, if the pac
option is specified. The phase and
amplitude of the signal is estimated using wavelets, with center
frequencies of fc
and fc2
for phase and amplitude components
respectively, where fc
is expected to be of lower frequency than
fc2
. Multiple frequencies for either the phase (fc
) or amplitude
(fc2
) can be specified with a commadelimited list, e.g. fc=1,2,4
.
Alternatively, a grid of frequencies cab be specified by using
fcrange=1,20
instead of fc
; here, it is also required to specify
num
, e.g. num=20
, which generates 20 fc
values evenly spaced on
a log scale, between 1Hz and 20Hz. The linear
option will force the grid
to be uniform on a linear, not log, scale instead. When considering phaseamplitude
coupling, only pairs where the phase frequency (fc
) is lower than
the amplitude frequency (fc2
) are retained for analysis.
As well as the wavelet center frequency, one can specify
the bandwidth of the wavelet, via the wavelet timedomain full width
at half maximum (FWHM), following this
convention.
See the CWTDESIGN
command for more
details, i.e. to better understand the implication of setting a
particular FWHM for the frequency domain properties of the wavelet.
In general, smaller FWHM values (in the timedomain) correspond to
broader FWHM in the frequency domain (i.e. a wider range of
frequencies above and below the center frequency will be captured).
Default FWHM values
If no FWHM values are specified, Luna will use a default value based on the frequency of the wavelet. For this, we seed on the default values selected in the Cox & Fell manuscript described in this vignette. In that manuscript, they used wavelet center frequencies of 0.5 to 30 Hz, evenly spaced on a log space, and paired these with 35 FWHM values evenlyspaced on a log scale from 5 to 0.25. For an arbitrary frequency, Luna will therefore use the following relation to generate a reasonable default value for the FWHM, which maintains equal frequency domain FWHM values on a logscale:
exp(0.7316762 * ln(F) + 1.1022791 )
Parameters
Core parameters
Parameter  Example  Description 

sig 
sig=C3,C4 
Optionally specify channels (default is to include all) 
fc 
fc=11,15 
Specify wavelet center frequency/frequencies (phase) 
fwhm 
fwhm=1,1 
Wavelet FWHM value(s) (phase) 
fc2 
fc2=11,15 
For PAC: as fc for amplitude 
fwhm2 
fwhm2=1,1 
For PAC: as fwhm for amplitude 
nreps 
nreps=1000 
Number of randomized, surrogate time series generated 
pac 
pac 
Estimate withinchannel phaseamplitude coupling metrics 
xch 
xch 
Estimate betweenchannel connectivity metrics 
Secondary parameters (primarily for specifying grids of frequency/FWHMs to test)
Parameter  Example  Description 

fcrange 
fcrange=1,20 
Specify range of center frequencies (requires num ) (phase) 
fwhmrange 
fwhmrange=5,0.25 
Specify range of FWHM values (requires num ) (phase) 
num 
num=20 
Number of intervals (evenly spaced on log scale) for fcrange or fwhmrange (phase) 
fcrange2 
fcrange2=1,20 
For PAC: as fcrange2 for amplitude 
fwhmrange2 
fwhmrange2=5,0.25 
For PAC: as fwhmrange2 for amplitude 
num2 
num2=20 
For PAC: as num for amplitude 
length 
length=10 
Wavelet duration (in seconds) 
linear 
linear 
For ranges (e.g. fcrange ) evenly space on a linear scale 
noepochoutput 
noepochoutput 
Suppress all epochlevel output 
That is, fcrange
and num
would be used instead of fc
, etc.
Output
Wholesignal level metrics (strata: CH1
x CH2
x F1
x F2
)
Variable  Description 

CFC 
0/1 for whether this metric is a crossfrequency contrast (i.e. F1 != F2 ) 
XCH 
0/1 for whether this metric is a crosschannel contrast (i.e. CH1 != CH2 ) 
dPAC 
Debiased phaseamplitude coupling metric 
dPAC_Z 
Empirical Zscore for dPAC based on timeshuffling 
wPLI 
Weighted phaselag index connectivity metric 
wPLI_Z 
Empirical Zscore for wPLI based on timeshuffling 
Epochlevel metrics (strata: E
x CH1
x CH2
x F1
x F2
)
Variable  Description 

dPAC 
Debiased phaseamplitude coupling metric 
dPAC_Z 
Empirical Zscore for dPAC based on timeshuffling 
wPLI 
Weighted phaselag index connectivity metric 
wPLI_Z 
Empirical Zscore for wPLI based on timeshuffling 
Example
See this vignette for an application of both dPAC and wPLI metrics.
Here, using two EEG channels from the second individual in the tutorial dataset, we extract only N2 sleep at consider wPLI between the two channels at a grid of 50 frequencies, linearly spaced between 1 and 25 Hz:
luna s.lst 2 o out.db
s 'MASK ifnot=NREM2 & RE &
CC sig=EEG,EEG(sec) fcrange=1,25 num=50 nreps=200 linear xch'
Extracting the output as follows:
destrat out.db +CC r F1 F2 CH1 CH2
We can plot the wPLI (left) and correspondong Zscore (the empirical value derived from randomization, right):
That is, we see significant connectivity between these two central channels in the spindle/sigma frequency band during N2 sleep.
Scaling the number of tests
Although this is designed to run reasonably efficiently, if you have lots of channels and/or lots of
frequencies and want to look at crossfrequency coupling between all pairs, then run time can get slow (...very slow...)
especially if you are using randomization to estimate Z scores (nreps
) and/or have lots of individuals in the datasets.
Therefore, try to focus these analyses and start small...
PSI
Calculates the phase slope index across channels
Estimates the phase slop index between pairs of channels, and provides a singlechannel summary of net PSI (i.e. net sender versus recipient).
The command requires that channels have similar sampling rates.
Parameters
To set the frequency or frequencies at which to calculate PSI either:

specify
f
andw
 one or more frequencies, for bandsfw/2
tof+w/2

or specify
flwr
,fupr
only  for a band of lower to upper (or an equal range of commadelimited values for each to give a range of bands, e.g.flwr=1,4,8
andfupr=4,8,12
implies three bands: 14 Hz, 48 Hz and 812 Hz) 
or specify
flwr
,fupr
,w
andr
 for a range of bands, fromf
equallingflwr
up tofupr
, in steps ofr
Hz; for each frequency, the window ofw
spanning the center frequency is constructed
Parameter  Example  Description 

sig 
sig=C3,C4 
An optional parameter, to specify which channels to calculate PSI 
epoch 
Report epochlevel (e.g. 30second epoch) output  
flwr 
flwr=3 
Consider frequencies from flwr to fupr (as a band, or in steps of r ) 
f 
f=10,12,15 
One or more frequencies to test (Hz) 
w 
w=5 
Window around each center frequency (Hz) in which to calculate PSI 
Secondary parameters are:
Parameter  Example  Description 

eplen 
eplen=5 
PSI subepoch length (default 4 seconds) 
seglen 
seglen=2.5 
Segment length (with 50% overlap) within each subepoch 
cachemetrics 
cachemetrics=c1 
Cache PSI, e.g. for use with PSC 
Cache options
Parameter  Example  Description 

cachemetrics 
cachemetrics=c1 
Cache net and pairwsie PSC (e.g. for PSC ) 
Output
Analysis parameter output (strata: F
)
Variable  Description 

F1 
Lower frequency 
F2 
Higher frequency 
NF 
Number of frequency bins 
Channellevel output (strata: CH
)
Variable  Description 

PSI 
Net PSI (standardized) for this channel 
Channel pair output (strata: CH1
x CH2
)
Variable  Description 

PSI 
Standardized PSI for this channel pair 
PSI_RAW 
Raw PSI 
STD 
Standard deviation of PSI 
Channellevel output (option: epoch
, strata: E
x CH
)
Variable  Description 

PSI 
Net PSI (standardized) for this channel 
Channel pair output (olption: epoch
, strata: E
x CH1
x CH2
)
Variable  Description 

PSI 
Standardized PSI for this channel pair 
PSI_RAW 
Raw PSI 
STD 
Standard deviation of PSI 
Example
to be added
MI
Calculates pairwise mutual information metrics across channels
Estimates mutual information, a measure of dependence between two signals, based on methods described in Analyzing Neural Time Series Data by MX Cohen. Total correlation and dual total correlation are two normalized variants of the mutual information statistic: MI / min[ H(X), H(Y) ] and MI/H(X,Y) respectively, where MI is mutual information, H(X) and H(Y) are the marginal entropies, and H(X,Y) is the joint entropy.
Signals are first coarsegrained: the number of bins is determined by one of three rules: FreedmanDiaconis (default), Scott or Sturges rule, as described in Cohen.
Parameters
Parameter  Example  Description 

sig 
sig=C3,C4,F3,F4 
Optionally specify channels (default is to include all) 
epoch 
epoch 
Calculate and report MI and other measures per epoch 
scott 
scott 
Use Scott's rule to determine bin number 
sturges 
sturges 
Use Sturges' rule to determine bin number 
permute 
permute=1000 
Estimate empirical significance via permutation, e.g. with 1000 null replicates 
Alert
Permutation and epochlevel analyses with the MI
command can be relatively slow .
Output
Output for the whole signal (strata: CH1
x CH2
)
Variable  Description 

MI 
Mutual information 
TOTCORR 
Total correlation 
DTOTCORR 
Dual total correlation 
JINF 
Joint entropy 
INFA 
Marginal entropy of first signal 
INFB 
Marginal entropy of second signal 
NBINS 
Number of bins 
Output perepoch, with epoch
(option: epoch
, strata: CH1
x CH2
)
Variable  Description 

MI 
Mutual information 
TOTCORR 
Total correlation 
DTOTCORR 
Dual total correlation 
JINF 
Joint entropy 
INFA 
Marginal entropy of first signal 
INFB 
Marginal entropy of second signal 
Output for permutation test (option: permute
)
Variable  Description 

EMP 
Empirical pvalue for mutual information statistic 
Z 
Zstatistic 
TSYNC
Crosscorrelation and phase delay
Estimate the crosscorrelation between two signals, within a window of W seconds, and report the estimated phase delay (in seconds) based on the maximal crosscorrelation in that time window.
Options
Parameter  Example  Description 

sig 
sig=C3,C4,F3,F4 
Optionally specify channels (default is to include all) 
w 
w=0.5 
Required time window (seconds) 
verbose 
Verbose output  
epoch 
Epochlevel output 
Output
Channelpair output (strata: CH1
x CH2
)
Variable  Description 

S 
Phase delay based on crosscorrelation 
Channelpair output (option: verbose
, strata: D
x CH1
x CH2
)
Variable  Description 

XCORR 
Estimated crosscorrelation for this delay 
Epochlevel channelpair output (option: epoch
, strata: CH1
x CH2
)
Variable  Description 

S 
Phase delay based on crosscorrelation 
Example
to be added