[Eeglablist] Rank deficiency and ASR+ICA

Fri Sep 27 12:53:06 PDT 2024

Hi Cedric,

> But I haven't found any perfect solution personally.

The perfect solution, if you want, is to process the EEG data in two
different ways.

1. High-pass filter at 1-2 Hz + ICA; This preprocessing is optimum for
ICA + time-frequency decomposition above the cutoff frequency of the HPF,
but not for time-domain averaged ERP because some people say
this 'distorts' the waveforms.
2. High-pass filter at < 0.5 Hz + ICA; This one is optimum for time-domain
averaged ERP but not so for ICA + time-frequency analysis.

Who said we can't do it, so why not?
I think using the above approach is more straightforward than worrying
about whether to use the HPF trick for ICA.

An important confirmation is to compare the ICA results between 1 and 2
(before that, repeat ICA several times on the same data to see how
reproducible the results are!) If they are practically the same, there is
no need to wonder about the choice.

Makoto

On Fri, Sep 27, 2024 at 1:50 PM Cedric Cannard via eeglablist <
eeglablist at sccn.ucsd.edu> wrote:

> Hi Francisco,
>
> I agree with Makoto it is not ideal solution.But I haven't found any
> perfect solution personally.
>
> If you do decide to go with the initial approach, here are the general
> steps.
>
> Yes, always get the best signal quality possible before ICA (although
> preserve eye movements as much as possible so they are extracted by ICA and
> not ASR, which removes relevant signals with the eye artifacts).
>
> A good way for this approach is to make a copy of your dataset before
> filtering. Filter copy data >1Hz, remove bad channels, run ASR, interpolate
> bad channels (after ASR to avoid data rank issues as there is no input like
> for ICA, to my knowledge), run ICA (input effective data rank). Store bad
> channels, bad segment latencies, ICA latencies during this process.
>
> Go back to the original dataset (<1 hz filter), remove bad channels,
> interpolate, remove bad segments, and transfer ICA weights (based on stored
> info). Run EEG = eeg_checkset(EEG); to make sure everything is good. Run
> ICLabel to flag and remove bad components.
>
> Hope this helps and that I didn't miss something.
>
> Cedric
>
>
>
>
> On Thursday, September 26th, 2024 at 2:59 PM, Cedric Cannard <
> ccannard at protonmail.com> wrote:
>
> > Hi Francisco,
> >
> > Yes, always get the best signal quality possible before ICA (although
> preserve eye movements as much as possible so they are extracted by ICA and
> not ASR, which removes relevant signals with the eye artifacts).
> >
> > A good way for this approach is to make a copy of your dataset before
> filtering. Filter copy data >1Hz, removebad channels, run ASR, interpolate
> bad channels (after ASR to avoid data rank issues as there is no input like
> for ICA, to my knowledge), run ICA (input effective data rank). Store bad
> channels, bad segment latencies, ICA latencies during this process.
> >
> >
> > Go back to original dataset (<1 hz filter), remove bad channels,
> interpolate, remove bad segments, and transfer ICA weights (based on stored
> info). Run EEG = eeg_checkset(EEG); to make sure everything is good. Run
> ICLabel to flag and remove bad components.
> >
> > Hope this helps and that I didn't miss something.
> >
> > Cedric
> >
> >
> >
> >
> >
> >
> > On Monday, September 23rd, 2024 at 10:13 AM, Francisco Zambrano via
> eeglablist eeglablist at sccn.ucsd.edu wrote:
> >
> > > Hello again,
> > >
> > > Thanks a lot for your insights, the results from
> > >
> > > Cedric
> > >
> > > 's reply gave me an
> > > answer that sounds more in line with what I've read.
> > >
> > > I've come with another question regarding ICA. If I'm trying to use
> ICA's
> > > trick of applying the >1 Hz data's unmixing matrix to the <1 Hz data,
> but
> > >
> > > I don't know if It's better to do the whole processing
> pipeline(removing
> > > bad channels, interpolating, removing bad trials,...) first and then
> > > applying it back to the <1Hz data that is also been processed with the
> same
> > > steps; or if I'm supposed to do it in another way. Could you please
> help me
> > > with this?
> > >
> > > Thanks in advance.
> > > Francisco.
> > >
> > > El mié, 7 ago 2024 a la(s) 3:37 p.m., Makoto Miyakoshi via eeglablist (
> > > eeglablist at sccn.ucsd.edu) escribió:
> > >
> > > > Hi Pancho,
> > > >
> > > > > rank(EEG.data(:,:)) over the original data It says the rank is 11
> > > >
> > > > Try rank(double(EEG.data(:,:)) and you'll see more realistic values
> at
> > > > least. It is counterintuitive if you see it naively that single or
> double
> > > > makes such a difference here.
> > > >
> > > > For the complete answer, see Cedric's reply.
> > > >
> > > > Makoto
> > > >
> > > > On Mon, Aug 5, 2024 at 9:40 PM Pancho Zambrano via eeglablist <
> > > > eeglablist at sccn.ucsd.edu> wrote:
> > > >
> > > > > Hi!,
> > > > >
> > > > > I'm currently working on a script for my thesis on an EEG channel
> > > > > reconstruction protocol for Universidad de Concepción, and I've
> found
> > > > > myself in some trouble regarding rank deficiency. The script it's
> about
> > > > > quantification of the positive yzing impact of ASR algorithm from
> > > > > cleanrawdata() on the data by analyzing the increase of the amount
> of
> > > > > dipolar sources present before and after using this default toolbox
> > > > > function. This is accomplished by comparing the residual variance
> from
> > > > > the
> > > > > dipfit model of the clean and original data. The problem is I use
> the
> > > > > matcorr function to find the best match between the original and
> clean
> > > > > ICA
> > > > > mixing matrices, but they have different number of columns: the
> first one
> > > > > is 64x64 and the clean one is 64x55. Since I've re-referenced the
> data
> > > > > after ASR, I thought this was the problem but I also noticed when
> I use
> > > > > rank(EEG.data(:,:)) over the original data It says the rank is 11,
> so I'm
> > > > > lost here, please help. I got this method from the article:
> > > > >
> > > > > Evaluation of Artifact Subspace Reconstruction for Automatic
> Artifact
> > > > > Components Removal in Multi-Channel EEG Recordings
> > > > >
> > > > > from:
> > > >
> > > >
> https://urldefense.com/v3/__https://ieeexplore.ieee.org/abstract/document/8768041?casa_token=IsR0oujDOEAAAAAA:T1NPgCRqT15sUYjG8KLtR8a9yW6PqQP7iaKzg5057R93Oef1yCmxDKt6ZtLM4IZk2t_xF4qWJ0I__;!!Mih3wA!ClbUhbvVFcY2DhNRP-wwMym8wBcZTXz2jeCl8yebysZiEmG01Z-cZm3OgV4T4R1GXwiDsWNEwyIFZcOglDZZkbBiLjHP0sFV$
> > > >
> > > > > This is my code:
> > > > > % Initialize EEGLAB
> > > > > eeglab;
> > > > >
> > > > > % Load the dataset
> > > > > EEG = pop_loadset('filename', 'C110_filt_butt_zap2.set');
> > > > >
> > > > > % Perform ICA on the original data
> > > > > EEG = pop_runica(EEG, 'extended', 1);
> > > > >
> > > > > % Save the original ICA results
> > > > > ica_weights_orig = EEG.icaweights;
> > > > > ica_sphere_orig = EEG.icasphere;
> > > > > ica_act_orig = EEG.icaact;
> > > > > ica_mix_matrix_orig = pinv(ica_weights_orig * ica_sphere_orig);
> > > > >
> > > > > % ASR dataset
> > > > > EEG_clean = pop_loadset('filename', 'prueba_1hz_ica.set');
> > > > >
> > > > > % % Perform ICA on the ASR cleaned data
> > > > > % EEG_clean = pop_runica(EEG_clean, 'extended', 1);
> > > > >
> > > > > % Save the ASR ICA results
> > > > > ica_weights_clean = EEG_clean.icaweights;
> > > > > ica_sphere_clean = EEG_clean.icasphere;
> > > > > ica_act_clean = EEG_clean.icaact;
> > > > > ica_mix_matrix_clean = pinv(ica_weights_clean * ica_sphere_clean);
> > > > >
> > > > > % Calculate the correlations between the components
> > > > > % n_components_orig = size(ica_mix_matrix_orig, 2);
> > > > > % n_components_clean = size(ica_mix_matrix_clean, 2);
> > > > > %
> > > > > % corr_matrix = zeros(n_components_orig, n_components_clean);
> > > > > %
> > > > > % for i = 1:n_components_orig
> > > > > % for j = 1:n_components_clean
> > > > > % corr_matrix(i, j) = corr(ica_mix_matrix_orig(:, i),
> > > > > ica_mix_matrix_clean(:, j));
> > > > > % end
> > > > > % end
> > > > >
> > > > > % Use the matcorr function to find the best matching
> > > > > rmmean = true; % Adjust as needed (remove mean or not)
> > > > > method = 0; % Correlation method
> > > > > weighting = []; % No specific weighting
> > > > >
> > > > > [corr, indx, indy, corrs] = matcorr(ica_mix_matrix_orig,
> > > > > ica_mix_matrix_clean, rmmean, method, weighting);
> > > > >
> > > > > % Display the matched correlations
> > > > > best_matches = corrs(sub2ind(size(corrs), indx, indy));
> > > > >
> > > > > disp('Correlations between the best matched components:');
> > > > > disp(best_matches);
> > > > >
> > > > > % Apply the spatial filter obtained from the ICA of the original
> data to
> > > > > the cleaned data
> > > > > W_orig = ica_weights_orig * ica_sphere_orig;
> > > > > Y_clean = W_orig * reshape(EEG_clean.data, size(EEG_clean.data,
> 1), []);
> > > > >
> > > > > % Calculate the activities of the ICs
> > > > > Y_clean = reshape(Y_clean, size(EEG_clean.icaact));
> > > > >
> > > > > % Calculate the mean power reduction for the IC activities
> > > > > Y_orig = W_orig * reshape(EEG.data, size(EEG.data, 1), []);
> > > > > Y_orig = reshape(Y_orig, size(EEG.icaact));
> > > > >
> > > > > power_reduction = mean(var(Y_orig, 0, 2)) - mean(var(Y_clean, 0,
> 2));
> > > > >
> > > > > disp('Mean power reduction:');
> > > > > disp(power_reduction);
> > > > >
> > > > > % Classify the ICs using iclabel
> > > > > EEG = iclabel(EEG, 'default');
> > > > > EEG_clean = iclabel(EEG_clean, 'default');
> > > > >
> > > > > % Get the IC classification labels
> > > > > ic_labels_orig = EEG.etc.ic_classification.ICLabel.classifications;
> > > > > ic_labels_clean =
> > > > > EEG_clean.etc.ic_classification.ICLabel.classifications;
> > > > >
> > > > > % Set up DIPFIT for the original data
> > > > > EEG = pop_dipfit_settings(EEG, 'hdmfile',
> > > > > 'standard_BEM/standard_vol.mat',
> > > > > ...
> > > > > 'coordformat', 'MNI', 'mrifile', 'standard_mri.mat', ...
> > > > > 'chanfile', 'standard_1005.elc', 'coord_transform', [0 0 0 0 0 0 1
> 1
> > > > > 1], ...
> > > > > 'chansel', 1:size(EEG.data, 1));
> > > > >
> > > > > % Fit dipoles to the original ICA components
> > > > > EEG = pop_multifit(EEG, 1:size(EEG.icaweights, 1), 'threshold',
> 100);
> > > > >
> > > > > % Get the residual variance of the fitted dipoles in the original
> data
> > > > > rv_orig = [EEG.dipfit.model.rv];
> > > > >
> > > > > % Set up DIPFIT for the ASR cleaned data
> > > > > EEG_clean = pop_dipfit_settings(EEG_clean, 'hdmfile',
> > > > > 'standard_BEM/standard_vol.mat', ...
> > > > > 'coordformat', 'MNI', 'mrifile', 'standard_mri.mat', ...
> > > > > 'chanfile', 'standard_1005.elc', 'coord_transform', [0 0 0 0 0 0 1
> 1
> > > > > 1], ...
> > > > > 'chansel', 1:size(EEG_clean.data, 1));
> > > > >
> > > > > % Fit dipoles to the ASR cleaned ICA components
> > > > > EEG_clean = pop_multifit(EEG_clean, 1:size(EEG_clean.icaweights,
> 1),
> > > > > 'threshold', 100);
> > > > >
> > > > > % Get the residual variance of the fitted dipoles in the ASR
> cleaned data
> > > > > rv_clean = [EEG_clean.dipfit.model.rv];
> > > > >
> > > > > % Identify dipolar ICs (residual variance < 5%)
> > > > > dipolar_ICs_orig = find(rv_orig < 0.05);
> > > > > dipolar_ICs_clean = find(rv_clean < 0.05);
> > > > >
> > > > > % Compare the number of dipolar ICs before and after ASR
> > > > > num_dipolar_ICs_orig = length(dipolar_ICs_orig);
> > > > > num_dipolar_ICs_clean = length(dipolar_ICs_clean);
> > > > >
> > > > > % Display dipolar IC results
> > > > > fprintf('Number of dipolar ICs in original data: %d\n',
> > > > > num_dipolar_ICs_orig);
> > > > > fprintf('Number of dipolar ICs in ASR cleaned data: %d\n',
> > > > > num_dipolar_ICs_clean);
> > > > >
> > > > > % Save all results
> > > > > save('ica_evaluation_results.mat', 'ica_weights_orig',
> > > > > 'ica_weights_clean',
> > > > > 'ica_sphere_orig', 'ica_sphere_clean', 'ica_act_orig',
> 'ica_act_clean',
> > > > > 'indx', 'indy', 'best_matches', 'Y_clean', 'power_reduction',
> > > > > 'ic_labels_orig', 'ic_labels_clean', 'rv_orig', 'rv_clean',
> > > > > 'dipolar_ICs_orig', 'dipolar_ICs_clean');
> > > > >
> > > > > Please, help.
> > > > >
> > > > > Francisco Alonso Zambrano Salamanca
> > > > > Biomedical Engineering
> > > > > Universidad de Concepción
> > > > > Chile
> > > > > _______________________________________________
> > > > > Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> > > > > To unsubscribe, send an empty email to
> > > > > eeglablist-unsubscribe at sccn.ucsd.edu
> > > > > For digest mode, send an email with the subject "set digest mime"
> to
> > > > > eeglablist-request at sccn.ucsd.edu
> > > > > _______________________________________________
> > > > > Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> > > > > To unsubscribe, send an empty email to
> > > > > eeglablist-unsubscribe at sccn.ucsd.edu
> > > > > For digest mode, send an email with the subject "set digest mime"
> to
> > > > > eeglablist-request at sccn.ucsd.edu
> > >
> > > _______________________________________________
> > > Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> > > To unsubscribe, send an empty email to
> eeglablist-unsubscribe at sccn.ucsd.edu
> > > For digest mode, send an email with the subject "set digest mime" to
> eeglablist-request at sccn.ucsd.edu
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to
> eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to
> eeglablist-request at sccn.ucsd.edu