[Eeglablist] Rank deficiency and ASR+ICA

[Makoto Miyakoshi]

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