spatio_temporal_cluster_test with NaNs in data

💬 Support & Discussions
,
1

Hey all,
new to MNE, and trying to use MNE’s spatio_temporal_cluster_test for some fNIRS data that has been separately processed using personal scripts (i.e., not using MNE etc).

As part of the pre-processing stage, artifacts are marked a NaNs, and I’m not sure what MNE/spatio_temporal_cluster_test is doing with the NaNs.

Artificially setting NaN to 0 across the data results in meaningful-looking results, but leaving NaNs in returns 0 clusters.

Is there any way of getting around having NaNs in the data?

For reference, I’m using MNE v 0.24.1 on Mac OS 12.3.1

Thanks in advance!
Mo

2

Hello @M055 and welcome to the forum!

I’m not quite sure where exactly you have NaNs in your data. Could you please give an example and maybe share a snippet of code that demonstrates what you’re doing?

Thanks!
Richard

3

Hello @richard
Sorry for the delay - I have NaNs because I have already-processed fNIRS [Oxy] data from an older dataset which I got as a 951(samples) X 24(channels) X 31(Subjects) from two conditions, and this data contains NaNs.

I’m basically trying to force previous non-MNE-data that contains NaNs into MNE’s stats.permutation_cluster_test function.

I did a little digging and I see that the nan_policy=‘omit’ is not yet implemented for the F-test in scipy, so it looks like my options might be to systematically remove NaNs from this data, or re-implement and pass an F-test that is ok with NaNs in the data.

I paste below a minimal sequence.

Thanks
Mo

import scipy.io
import numpy as np
from mne import stats

# Adjacency matrix from another file:
import NIRS_adj_mat
adj_matrix = NIRS_adj_mat.get_NIRS_adj_mat()

# Load the data from the obtained mat file
matdata = scipy.io.loadmat('./matlab_data/avg_keep.mat')

# Recover data for the 'R' and 'N' conditions
allOxyR = np.array(matdata['avg'][0])
allOxyN = np.array(matdata['avg'][1])

allOxyR = np.moveaxis(allOxyR,-1,0)
allOxyN = np.moveaxis(allOxyN,-1,0)

# Shape of arrays is 31(subjects) X 951(samples) X 24(channels)
allOxyR.shape, allOxyN.shape
>> ((31, 951, 24), (31, 951, 24))

# Make a copy where NaNs are replaced by 0
allOxyR0=allOxyR.copy()
allOxyN0=allOxyN.copy()
allOxyR0[np.where(np.isnan(allOxyR0))]=0
allOxyN0[np.where(np.isnan(allOxyN0))]=0

# Create an 'X' matrix with original data containing NaN and one with NaNs replaced by 0
X = [allOxyR, allOxyN]
X0 = [allOxyR0, allOxyN0]

# Run MNE clustering
cluster_stats = stats.permutation_cluster_test(X,adjacency=adj_matrix)
>> Using a threshold of 4.001191
>> stat_fun(H1): min=nan max=nan
>> Running initial clustering
>> Found 0 clusters

cluster_stats = stats.permutation_cluster_test(X0,adjacency=adj_matrix)
>> Using a threshold of 4.001191
>> stat_fun(H1): min=nan max=nan
>> Running initial clustering
>> Found 42 clusters
>> Permuting 1023 times...
4

Hi @M055 ,
I am trying to perform spatio temporal cluster tests with my fNIRS data and I encounter the problem of the adjacency, since nirs channels are not implemented yet for the ch_find_adjacency() function.
I see in your code that you used a personnalized function called NIRS_adj_mat.
Could tell me more about it? How did you obtain an adjacency matrix?
I am using a NIRx system and I imported my info object from a raw NIRx folder.
Thanks so much !
Jérémie

5

Hi Jeremie,
I just manually made a list of all adjacent channel pairs of the 12LH+12RH, used that to fill in 1s and 0s in a 24x24 numpy matrix, and converted it to a boolean(np.bool8) , sparse(scipy.sparse) matrix and saved it. Please see example code below.
Best,
Mo

Edit NB: Not sure why channel numbers on the adjacency matrix are not 0=indexed (we moved away from MNE a while ago)

def get_NIRS_adj_mat():
    '''
    return adjacency matrix
    
    '''
    import numpy as np
    import scipy.sparse as sps
    connected_channels = [1, 3, 
                         1, 4,
                         2, 4,
                         ...
                         ...
                         ]
    connected_channels = np.array(connected_channels).reshape( int(len(connected_channels)/2),2)
    ch_adj_mat = np.zeros((24,24))
    for x,y in zip(connected_channels[:,0],connected_channels[:,1]):
        ch_adj_mat[x,y]=1
    ch_adj_mat = np.bool8(ch_adj_mat)
    ch_adj_mat_sparse = sps.csr.csr_matrix(ch_adj_mat)

    return ch_adj_mat_sparse
6

Hi Mo,
Thank you so much for your answer, in the meantime, I did the same as you but by taking as an example a fieldtrip-like neighbour layout and it worked !
Have a good one :slight_smile:
Jérémie

1 Like