ANOVA Group analysis

💬 Support & Discussions
1

Dear MNE team,

I would like to perform a 2x2 ANOVA with one between-group variable (Patient vs Control) and one within-group variable (Child vs Adult). I am encountering an issue with assigning my conditions to my data. Here is the code for assignation:

factor_levels = [2, 2]
effects = ['A', 'B', 'A*B']

or ri, r in enumerate(config.group_anova):
    for di, d in enumerate(config.age):

        if ri == 0 and di == 0:
            cidx = 0
        else:
            cidx = cidx + 1

        for si, sid in enumerate(sid_list):

            try:
                print(f'Processing: {sid}')
                evo_in = f'{config.ievo_dir}{sid}_age_{d}_{r}-ave.fif'
                tmp = mne.read_evokeds(evo_in)[0].apply_baseline(baseline=bsln).crop(tmin=stat_tmin, tmax=stat_tmax)
                tmp.comment = f'{sid} {cnames[cidx]}'

                data[cidx][si] = np.transpose(tmp.data)

In my data matrix, I have values of 0 because both factors are being considered as within-group variables (when one is actually a between-group factor). Could you please help me resolve this issue with the assignment of within- and between-group variables? Thank you in advance.

data

    ...,
     [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,
       0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
     [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,
       0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
     [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,
       0.00000000e+00,  0.00000000e+00,  0.00000000e+00]],
    ...,
    [[ 2.93491243e-07,  8.77661356e-08, -1.49599865e-07, ...,
      -9.20545059e-08, -9.00147917e-07, -5.93777352e-07],
     [ 2.24999015e-07,  1.28827063e-08, -1.51916348e-07, ...,
       3.63874059e-08, -8.80130666e-07, -6.97680029e-07],
     [ 9.78560195e-08, -8.65436198e-08, -1.05204335e-07, ...,
       1.74845550e-07, -8.05547142e-07, -7.05786174e-07],
     ...................................
  • MNE version: 0.17.1
  • operating system: Windows 11
2

what exactly is your sid_list?

From the two for loops, it seems that you want to only read the subjects that are either children/adults and/or patients/controls. However, you use the same sid_list every time.

Furthermore: how do you initialize data? Is it a python list or a numpy array?

3

Hello Marijn, thank you for your response. Here is the initialization of my variables (data is a numpy array)

patients_sids = [f'S{i}' for i in range(1, 10)]  # S1 to S11
controls_sids = [f'S20{i}' for i in range(1, 10)]  # S201 to S209
sid_list = patients_sids + controls_sids  # Combined list of subjects

cmap = plt.cm.Greys
cols = ('lightgrey', 'darkgrey', 'lightgrey', 'darkgrey')
ext = f'_age'

n = len(controls_sids)  # Number of controls (assuming same number of subjects in both groups)
stat_dir = config.stat_dir + 'Evo_Anova/'

# Naming convention
cnames = ['P_Child', 'P_Adult', 'C_Child', 'C_Adult']
stat_tmin, stat_tmax = 0.0, 0.8
bsln = (-.2, 0)

#_____________  Initiate report and include av data plot
report_out  = f'{config.stat_dir}/EmoAge_anova/EmoAge_anova_anova.html'
report      = Report(report_out, verbose=False)

factor_levels = [2, 2]
effects = ['A', 'B', 'A:B']
return_pvals = False


# ANOVA function
def stat_fun(*args):
    # get f-values and p-values
    return f_mway_rm(np.swapaxes(args, 1, 0),
                     factor_levels=factor_levels,
                     effects=effects,
                     return_pvals=return_pvals)[0]

#%% ______________  [3] Load gd av for vizualisation
evo_list = list()
com_list = list()`
4

I found another way to assign conditions that seems to work

# Concatenation 
event_id = {'1': 1, '2': 2, '3': 3, '4': 4, '5': 5, '6': 6}
combined_epochs = mne.concatenate_epochs(Patient_Neutral_list + Patient_Disgust_list + Patient_Happy_list + Control_Neutral_list + Control_Disgust_list + Control_Happy_list)
combined_epochs.equalize_event_counts(event_id)


# Data transformation
X = [combined_epochs[event_name].get_data() for event_name in event_id]
X = [np.transpose(x, (0, 2, 1)) for x in X]

# Adjency
adjacency, ch_names = find_ch_adjacency(combined_epochs.info, ch_type="eeg")

# ANOVA
n_subjects = len(sid_patients_list) + len(sid_control_list)
factor_levels = [2, 3]
effects = ['A', 'B', 'A:B'] #'A*B'
f_thresh = f_threshold_mway_rm(n_subjects, factor_levels, effects, pthresh)

# Cluster test
for fi, f in enumerate(f_thresh):
    T_obs, clusters, p_values, H0 = cluster_stats = spatio_temporal_cluster_test(
        X, adjacency=adjacency, threshold=f, n_permutations=n_permutations, tail=1
    )

    good_cluster_inds = np.where(p_values < p_accept)[0]
    print(f'Clusters trouvés: {good_cluster_inds}')