[Eeglablist] Rank deficiency and ASR+ICA

[Pancho Zambrano]

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