Creating TemporalDataset object from subjects with different channels

[jimzhang629]

Hi!

I'm trying to use this lovely toolbox for intracranial EEG analysis. However, my subjects have different numbers of channels, and different channel locations/names. I've already loaded my data into a nested dictionary that's structured in terms of subject -> condition -> Epochs object. However, when I try to concatenate across Epochs objects, I get a ValueError because different subjects have different numbers of channels. This may be a silly question but is there any way to have all of my subjects in the same TemporalDatasets object? I would like to have them in the same dataset so that I can group electrodes in the same ROI across subjects. Thank you in advance for any advice!

import mne
import numpy as np
from rsatoolbox.data import TemporalDataset

def create_temporal_dataset(data_dict):
    """
    Create a TemporalDataset from a dictionary of Epochs data.
    
    Parameters:
    data_dict (dict): Dictionary with the structure subject -> condition -> HG_ev1_rescaled.
    
    Returns:
    TemporalDataset: The rsatoolbox TemporalDataset object created from the dictionary data.
    """
    # Initialize lists to hold the data and labels
    data = []
    obs_descriptors = {'subject': [], 'condition': []}
    channel_descriptors = {}
    time_descriptors = {}

    # Extract the data
    for subject, conditions in data_dict.items():
        for condition, epochs_dict in conditions.items():
            epochs = epochs_dict['HG_ev1_rescaled']
            # Get the data from the Epochs object
            epochs_data = epochs.get_data()  # shape: (n_epochs, n_channels, n_times)
            
            # Append data and descriptors
            data.append(epochs_data)
            obs_descriptors['subject'].extend([subject] * epochs_data.shape[0])
            obs_descriptors['condition'].extend([condition] * epochs_data.shape[0])

    # Convert lists to numpy arrays
    data = np.concatenate(data, axis=0)  # shape: (total_epochs, n_channels, n_times)

    # Assuming the first Epochs object has the relevant info for channels and times
    sample_epochs = next(iter(next(iter(data_dict.values())).values()))['HG_ev1_rescaled']
    channel_descriptors['name'] = sample_epochs.ch_names
    time_descriptors['times'] = sample_epochs.times

    # Create TemporalDataset
    dataset = TemporalDataset(
        measurements=data,
        obs_descriptors=obs_descriptors,
        channel_descriptors=channel_descriptors,
        time_descriptors=time_descriptors
    )

    return dataset

dataset = create_temporal_dataset(subjects_mne_objects)
print(dataset)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[8], [line 52](vscode-notebook-cell:?execution_count=8&line=52)
     [49](vscode-notebook-cell:?execution_count=8&line=49)     return dataset
     [51](vscode-notebook-cell:?execution_count=8&line=51) # Example usage:
---> [52](vscode-notebook-cell:?execution_count=8&line=52) dataset = create_temporal_dataset(subjects_mne_objects)
     [53](vscode-notebook-cell:?execution_count=8&line=53) print(dataset)

Cell In[8], [line 34](vscode-notebook-cell:?execution_count=8&line=34)
     [31](vscode-notebook-cell:?execution_count=8&line=31)         obs_descriptors['condition'].extend([condition] * epochs_data.shape[0])
     [33](vscode-notebook-cell:?execution_count=8&line=33) # Convert lists to numpy arrays
---> [34](vscode-notebook-cell:?execution_count=8&line=34) data = np.concatenate(data, axis=0)  # shape: (total_epochs, n_channels, n_times)
     [36](vscode-notebook-cell:?execution_count=8&line=36) # Assuming the first Epochs object has the relevant info for channels and times
     [37](vscode-notebook-cell:?execution_count=8&line=37) sample_epochs = next(iter(next(iter(data_dict.values())).values()))['HG_ev1_rescaled']

ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 175 and the array at index 16 has size 174

[JasperVanDenBosch]

Hi! Great to hear that you're considering using this on intracranial data.

I think one solution would be to add NaNs? i.e. make a master array that has all channels, before creating the Dataset. For RDMs there is rsatoolbox.rdm.combine.from_partials which will do this for you, but I don't think we have something for Datasets yet.

[JasperVanDenBosch]

There's also merge_datasets but this assumes channels are the same across the partitions.

[jimzhang629]

Thanks! Yeah it seems like adding NaNs is the best solution.

I also want to add ROI and subject information for each channel, I suppose I can just add those as lists similar to how the channel names are added as a list in the demo_temporal.ipynb? And then just have new keys of 'ROI' and 'subject' in the channel_descriptors dictionary?

So it would be something like:

channel_names = ['CH1', 'CH2', 'CH3']
roi = ['v1', 'v1', 'v2']
subjects = ['sub1', 'sub1', 'sub1']

chn_des = {'channels': channel_names, 'roi': roi, 'subjects': subjects}

data = rsatoolbox.data.TemporalDataset(measurements, 
                              descriptors = des, 
                              obs_descriptors = obs_des, 
                              channel_descriptors = chn_des, 
                              time_descriptors = tim_des)

Or would a better approach to be have a separate TemporalDataset for each roi? This is because in the master channel array, for each channel, if I fill in trials with NaNs on which that channel doesn't actually have data, then I'll end up with a huge array. I'm just not sure if having separate TemporalDatasets would make doing the RSA harder.

[JasperVanDenBosch]

This looks good except that subject would be an observation descriptor

[jimzhang629]

Thank you for all of your help so far! I managed to make a TemporalDataset for each roi. However, my issue now is that when I try to use the calc_rdm_movie function, I get all NaNs for my dissimilarities. I suspect this is because I have NaNs on each trial now.

Here's exactly what I did in case more context is needed. First I prepare my data so I can convert it to a TemporalDataset object:

import numpy as np

def prepare_data_for_temporal_dataset(subjects_mne_objects, condition_names, rois, subjects, sig_electrodes_per_subject_roi):
    # returns dat
    dat = {}
    overall_electrode_mapping = []
    electrode_mapping_per_roi = {roi: [] for roi in rois}  # Reinitialize for each processing run
    print('subjects: ', subjects)
    for roi in rois:
        dat[roi] = {}  # make a dict for each roi
        dat[roi]['channel_names'] = []  # initialize a list to hold channel names
        dat[roi]['channel_rois'] = [] # initialize a list to hold what roi each channel is a part of
        dat[roi]['condition_names'] = {}  # initialize dict where keys are condition names and values are integer indices
        dat[roi]['cond_idx'] = np.array([], dtype=int)  # initialize an empty 1D array for condition indices for each trial
        dat[roi]['sub_idx'] = np.array([], dtype=int) # initialize an empty 1D array for subject for each trial
        dat[roi]['times'] = np.array([])  # initialize an empty 1D array for time points

        # Determine all unique channels across subjects for this ROI, maintaining order
        all_channels = []
        for sub in subjects:
            sig_electrodes = sig_electrodes_per_subject_roi[roi].get(sub, [])
            sub_channel_names = [sub + '-' + sig_electrode for sig_electrode in sig_electrodes]
            for chan in sub_channel_names:
                if chan not in all_channels:
                    all_channels.append(chan)
                    dat[roi]['channel_rois'].append(roi) # append roi for each channel
        dat[roi]['channel_names'] = all_channels
        num_channels = len(all_channels)
        print('num channels: ', num_channels)

        # Initialize the data array with the number of trials and total channels
        dat[roi]['data'] = np.empty((0, num_channels, 0))  # initialize an empty 3D array for trials x channels x time points

        for sub in subjects:
            total_trials = 0  # Initialize counter for total trials across all conditions
            sig_electrodes = sig_electrodes_per_subject_roi[roi].get(sub, [])
            sub_channel_names = [sub + '-' + sig_electrode for sig_electrode in sig_electrodes]
            if not sig_electrodes:
                continue

            cond_idx = 0  # the example uses indexing from 1, but let's start from 0 because python
            for condition_name in condition_names:
                print(f'Processing {sub} for {condition_name} in {roi}')
                epochs = subjects_mne_objects[sub][condition_name]['HG_ev1_rescaled'].copy().pick(sig_electrodes)
                dat[roi]['condition_names'][condition_name] = cond_idx

                epochs_data = epochs.get_data(copy=True)
                num_trials, num_sub_channels, num_timepoints = epochs_data.shape

                print(f'Number of trials for {sub} in {condition_name}: {num_trials}')
                total_trials += num_trials

                # Initialize data array time dimension if it is empty
                if dat[roi]['data'].shape[2] == 0:
                    dat[roi]['data'] = np.empty((0, num_channels, num_timepoints))
                    dat[roi]['times'] = epochs.times

                # Create an array filled with NaNs for the current subject's data
                sub_data = np.full((num_trials, num_channels, num_timepoints), np.nan)

                # Find the indices for the subject's channels in the total list of channels
                channel_indices = [all_channels.index(chan) for chan in sub_channel_names]
                print('sub: ', sub)
                print("channel indices: ", channel_indices)
                
                # Place the subject's data in the correct indices
                sub_data[:, channel_indices, :] = epochs_data

                # Concatenate the new data along the first axis (trials)
                dat[roi]['data'] = np.concatenate((dat[roi]['data'], sub_data), axis=0)

                # Extend the cond_idx array
                dat[roi]['cond_idx'] = np.concatenate((dat[roi]['cond_idx'], np.full(num_trials, cond_idx)))

                # extend the sub_idx array by this number of trials with this subject
                dat[roi]['sub_idx'] = np.concatenate((dat[roi]['sub_idx'], np.full(num_trials, sub)))
                print('last entry in cond_idx is: ', dat[roi]['cond_idx'][-1])
                print('cond_idx is: ', cond_idx)
                cond_idx += 1  # increment cond_idx

            print(f'Total number of trials for {sub} across all conditions: {total_trials}')

    return dat

# call the function
dat = prepare_data_for_temporal_dataset(subjects_mne_objects, condition_names, rois, subjects, sig_electrodes_per_subject_roi)

Next, I make the temporal datasets

# now make temporal datasets from prepared data
temporal_datasets = {}

for roi in rois:
    measurements = dat[roi]['data']
    des = {'roi': roi}
    obs_des = {'cond_idx': dat[roi]['cond_idx'], 'sub_idx': dat[roi]['sub_idx']}
    chn_des = {'channel_names': dat[roi]['channel_names'], 'channel_rois': dat[roi]['channel_rois']}
    tim_des = {'time': dat[roi]['times']}

    temporal_datasets[roi] = rsatoolbox.data.TemporalDataset(measurements, 
                                             descriptors = des, 
                                             obs_descriptors = obs_des, 
                                             channel_descriptors = chn_des,
                                             time_descriptors = tim_des)

And when I look at my data:
temporal_datasets['occ'].sort_by('cond_idx')
I see that my data has NaNs on each trial, as expected.

rsatoolbox.data.TemporalDataset(
measurements = 
[[[-0.00600946  0.04795925  0.08308495 ... -0.69147283 -0.62073958
   -0.51093721]
  [        nan         nan         nan ...         nan         nan
           nan]
  [        nan         nan         nan ...         nan         nan
           nan]
  ...
  [        nan         nan         nan ...         nan         nan
           nan]
  [        nan         nan         nan ...         nan         nan
           nan]
  [        nan         nan         nan ...         nan         nan
           nan]]

 [[        nan         nan         nan ...         nan         nan
           nan]
  [        nan         nan         nan ...         nan         nan
           nan]
  [        nan         nan         nan ...         nan         nan
           nan]
  ...
  [        nan         nan         nan ...         nan         nan
           nan]
...

But this doesn't play nice with the calc_rdm_movie function:

rdms_occ_data = calc_rdm_movie(temporal_datasets['occ'], method='correlation', descriptor='cond_idx')

Which outputs all nans for the dissimilarities:

dissimilarity_measure = 
correlation
dissimilarities = 
[[nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 ...
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]]

[jimzhang629]

I think this is because the np.einsum step in the calc_rdm_correlation function doesn't work with nans, and will propagate them.

[JasperVanDenBosch]

That's right, for a Dataset with NaN's like that you'd have to use rsatoolbox.rdm.calc_rdm_unbalanced, see also this demo

[jimzhang629]

It worked! Thanks so much for all of your help, Jasper!