[Eeglablist] Rank deficiency and ASR+ICA

[Cedric Cannard]

Dear Pancho,

I don't know about these other pieces but for getting the correct effective data rank, use this instead of MATLAB's rank function:

data_rank = sum(eig(cov(double(TMPEEG.data')))>1E-7);  % continuous data
% data_rank = sum(eig(cov(double(TMPEEG.data(:,:)'))) > 1E-7); % epoched data


Cedric



On Monday, August 5th, 2024 at 4:49 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