Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

noise_prec_from_unbalanced produces NaN values #389

Open
ndolf opened this issue Feb 23, 2024 · 4 comments
Open

noise_prec_from_unbalanced produces NaN values #389

ndolf opened this issue Feb 23, 2024 · 4 comments

Comments

@ndolf
Copy link

ndolf commented Feb 23, 2024

Hello,

I recently started using rsatoolbox in python in order to compute Mahalanobis distances (and eventually crossnobis + taking noise variance into account). I have an unbalanced design - the number of voxels/channels differs between conditions and I am not sure how the toolbox is dealing with these NaN values when computing the precision matrix.

I have 8 different conditions (or items) and 4 runs, resulting in 32 patterns (Each condition was repeated once in every run). The number of channels differ between conditions (or even within conditions across runs) - when a value is missing this is replaced by a NaN value

            measurements = measurements['activation_patterns']
            nItems = 8
            n_rep = measurements.shape[0]/nItems;
            nCond = measurements.shape[0]/n_rep;       
            nVox = measurements.shape[1];
            
            items = np.array(['stim_%02d' % x for x in np.arange(nCond)])
            items = np.repeat(items, n_rep)
            runs = np.tile(np.arange(n_rep), nItems)
            
            # now create a  dataset object
            des = {'Task': task, 'subj': sub, 'ROI': roi_counter+1}
            obs_des = {'items': items, 'runs': runs}
            chn_des = {'voxels': np.array(['voxel_' + str(x) for x in np.arange(1,nVox+1)])}           
         
            data=rsd.Dataset(measurements=measurements,
                                       descriptors=des,
                                       obs_descriptors=obs_des,
                                       channel_descriptors=chn_des)

I wanted to start simple and follow the demo to compute the precision matrix:
noise_prec_diag = rsatoolbox.data.noise.prec_from_measurements(data, obs_desc='items', method='diag')

as this gave me a noise precision matrix with NAN values I tried:

         noise_prec_diag = rsatoolbox.data.noise.prec_from_unbalanced(data, obs_desc='items', method='diag')

Again I get all NaN values in the noise precision matrix and I noticed when looking at data.measurements that running this function sets all the rows of any channel that had missing values to NaN (so now the channel has NaN values for all conditions and runs, while initially it only had NaN for certain rows).

I tried to figure it out by looking at the code of the toolbox but it is really hard to follow the order of computations.

I have two questions:

  • The precision matrix has 1 variance per voxel, so this means that the variance for each voxel is computed across all conditions and runs. Is this how I am supposed to do it?
  • Why do I get NaN in the output of the precision estimation? Does it not deal with missing data (even though it is the function for unbalanced designs)?

Irrespective of the issue with the NAN values I was also doubting whether 8 conditions x 4 runs is enough to estimate precision. I have residuals (1 per voxel) from the GLM that potentially could be used instead but I was not sure it is ok to just have 1 residual per voxel (even though I have 32 multivoxel patterns = 8 conditions x 4 runs).

Thank you in advance for any help!

Nina

@JasperVanDenBosch
Copy link
Contributor

Hi Nina,
Thanks for bringing this up. The Dataset should definitely not be changed by the noise function. I'm gonna look into that. Could you send me some data to test with? (Like with Dataset.save()). But it would make sense that the noise matrix does not have values for the channel which is Nan.
Why do you end up with only one residual value per voxel? Is the GLM on concatenated runs?

@ndolf
Copy link
Author

ndolf commented Mar 5, 2024

Hi Jasper,

thanks for writing back to me! In the attachment you can find one dataset object for one particular ROI (lateral occipital complex - 423 voxels).

sub01_roi2_localizer_Runs.zip

Indeed, I have all my runs in one single GLM - I use SPM. Below I include a example of my design matrix

example_GLM

Let me know if you have more questions.

Nina

@ndolf
Copy link
Author

ndolf commented Apr 16, 2024

Dear Jasper,

I am here at CNS listening to a talk of Niko Kriegeskorte about the rsatoolbox and this reminded me of my post. Did you have time to look at the issue I am having?

Thanks!

@JasperVanDenBosch
Copy link
Contributor

Hi Nina, unfortunately I haven't yet, but thanks for reminding!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants