<div>Thanks Jan for that detailed discussion.  I'd like to just focus on issues regarding the PLV, the topic of recent discussion.</div>
<div> </div>
<div>Could someone discuss the present status of the phase locking value? </div>
<div> </div>
<div>I recall that when Lachaux et al. (1999) introduced it (I think it was them) I didn't understand what is the advantage of adding unit magnitude phasors. I would have thought that neglecting the magnitude of the phasor would give a noisier measure of coherence.</div>


<div> </div>
<div>Another worry is that the PLV includes the real part of the phasor which is contaminated by the volume conduction from the forward model.  To avoid that contamination many people only look at the imaginary part (the aspect of PLV with 90 deg phase shift). The problem with the imaginary part is that one loses the information from the non-volume conduction real part.  </div>


<div> </div>
<div>One also needs to consider the potential contamination due to the reference. So PLV is a topic on which care must be taken. </div>
<div> </div>
<div>I'm interested in the current thinking on these issues.</div>
<div> </div>
<div>thanks</div>
<div>Stan</div><br><br>
<div class="gmail_quote">On Fri, Mar 25, 2011 at 12:08 PM, jan-mathijs schoffelen <span dir="ltr"><<a href="mailto:jan.schoffelen@donders.ru.nl">jan.schoffelen@donders.ru.nl</a>></span> wrote:<br>
<blockquote style="BORDER-LEFT: #ccc 1px solid; MARGIN: 0px 0px 0px 0.8ex; PADDING-LEFT: 1ex" class="gmail_quote">
<div style="WORD-WRAP: break-word">Hi Matt, Marco and others, 
<div><br></div>
<div>In the following I will address some issues which occurred earlier in this thread and try to also address the stuff that has come up since.</div>
<div><br></div>
<div>Here's my 2 cents:</div>
<div><br></div>
<div>As Matt and others already found out to their frustration, there is no straightforward way (in terms of how fieldtrip can handle the data) to carry out a statistical test (across a group of subjects) on data, which is the output of ft_connectivityanalysis. I will not talk about single subject analyses, but I sensed some confusion here and there, in that the 'statfun_indepsamplesZcoh' occasionally popped up. This is essentially indeed a function allowing to compute significance of coherence differences within a single subject.</div>


<div>Let's return to the statistics across a group case. The dirty hack I proposed a while ago in principle should work, but I could imagine that it doesn't work flawlessly, and it may need some fine-tuning. Having said this, in the past some functionality was built in to actually apply this dirty trick but only specifically so on coherence data, i.e. the data field called 'cohspctrm'. Internally in the ft_freqstatistics code, the field 'label' was created, in which each element consisted of a concatenation of 2 label strings, e.g. {'E1 - E2'}. The reason why other metrics are not supported, is that the aforementioned functionality was built in long before the existence of ft_connectivityanalysis, at the time we were mostly doing coherence analysis. As such, there is no principled reason why fieldtrip should not support similar functionality for other connectivity metrics, and it's certainly on the developers' to-do list (in this case:  my to-do list).</div>


<div>On a related note, once one manages to 'fool' freqstatistics (or if one would like to do a more traditional analysis with freqstatistics) and wants to apply clustering in the channel dimension, fieldtrip needs to know which channels are neighbours of each other. For this, the cfg structure needs to have a field called cfg.neighbours, and if this does not exist, fieldtrip tries to create one on the fly, based on the specification of the positions of the channels. This information is always stored (if known) in the data field called grad (for MEG), or elec (for EEG). Apparently, Marco's data does not contain this field, most likely causing the error he reported in an earlier message. This does not mean that one cannot do clustering in this case. It just means that one has to create the neighbourhood structure between the channels on his own. There's surely some documentation about this on the wiki, so I don't want to enter into details about this here.</div>


<div>If one now has scalp level estimates of connectivity between channel pairs, the clustering over channels may in principle be done (but not yet in fieldtrip as far as I am aware), since the channel dimension is now not just defined in 2D (as a projection of the channels onto a plane), but in 4-dimensional space. So in such a case one needs to specify an empty neighbourhood structure, i.e. cfg.neighbours{k}.label = data.label{k} and cfg.neighbours{k}.neighblabel = {} (where k is the k'th label in the (in this case) bivariate data, see above. Fieldtrip then still allows for clustering over frequencies and time.</div>


<div>Then it is of course the question as to how to generate the connectivity data in order to optimally fool ft_freqstatistics, and one level higher in the analysis pipeline, it relates to how ft_freqanalysis should behave. Let me start from the top, and use coherence as an example. Coherence essentially is the frequency domain version of the cross-correlation function and as such is constructed from the frequency domain representation of the time series. Analogous to the correlation which is the covariance between two signals, normalised with the variance of the two respective signals, the coherence is the cross-spectral density between two signals, normalised with the power of the two respective signals, at frequency X. The correlation can be loosely seen as a multiplication between two time domain signals, the cross-spectral density is obtained by (conjugate) multiplication between the frequency domain representation of those two signals. This essentially is related to the two different (relevant) outputs ft_freqanalysis will provide. Cfg.output = 'fourier' provides the fourier transforms of the time domain signals, and cfg.output='powandcsd' provides the cross-spectral density between a specified set of channel combinations, and you get the power spectra for free.</div>


<div>Both versions allow for the coherence to be computed but the memory requirements at this stage differs. For 'fourier' data the dimensionality of the fourierspctrm = 'rpttap_chan_freq_time' (by force the single replicates are kept in memory, because replicates may only be averaged once the fourierspectra are multiplied (yielding auto- and cross-spectra). For 'powandcsd' data the dimensionality of crsspctrm is either 'chan_freq_time', or 'rpt_chan_freq_time'. Averaging across replicates is now possible (when cfg.keeptrials = 'no' before calling ft_freqanalysis)  (since the multiplication is now performed), yet the 'chan' dimension now represents channel combinations, which in the case of cfg.channelcmb = {'all' 'all'} leads to (N-1)*N/2 channel combinations. This typically leads to a much higher memory usage, particularly because for some connectivity metrics (such as plv) one needs to process the single replicates before computing the plv, requiring cfg.keeptrials to be 'yes' prior to calling ft_freqanalysis. Therefore I would advocate the use of cfg.method = 'fourier', as long as the number of trials and tapers is not too big.</div>


<div>Next, ft_connectivityanalysis can be called, using the freq data with fourier. One can actually specify a cfg.channelcmb at this stage, which makes sense if not all combinations are needed, or if the requested connectivity metric is symmetric, such as coherence. This could be the trick Matt is waiting for, i.e. run ft_freqanalysis with output='fourier', and specify cfg.channelcmb = {'all' 'all'} before calling ft_connectivityanalysis. Alternatively, you could do a reshape on the plvspctrm, and keep only the interesting non-duplicated pairs, e.g. indx = reshape(1:length(conn.label).^2,[length(conn.label) length(conn.label]); indx = tril(indx,-1); plvspctrm = reshape(plvspctrm, [length(conn.label) length(conn.label) length(conn.freq) length(conn.time)]); plvspctrm = plvspctrm(indx(:),:,:); label = label(indx); label in this case is the processed labelcmb field.</div>


<div><br></div>
<div>I hope that this will be of any help.</div>
<div><br></div>
<div>Jan-Mathijs</div></div></blockquote></div>