lmer changes

examples/EEG_AD/figures/EEG_predictions.py

@@ -14,16 +14,14 @@
 # Set the p-value threshold
 p_value_th = 0.01
 
-def append_lmer_results(lmer_results, group, elec, p_value_th, data_lmer):
+def append_lmer_results(lmer_results, elec, p_value_th, data_lmer):
     '''
-    Create a list with the z-scores of the linear mixed model analysis for a given group and electrode.
+    Create a list with the z-ratios or t-ratios of the linear mixed model analysis for a given group and electrode.
 
     Parameters
     ----------
     lmer_results : dict
         Dictionary with the results of the linear mixed model analysis.
-    group : str
-        Group name.
     elec : int
         Electrode index.
     p_value_th : float
@@ -37,11 +35,15 @@ def append_lmer_results(lmer_results, group, elec, p_value_th, data_lmer):
         List with the z-scores of the linear mixed model analysis.
     '''
 
-    p_value = lmer_results[f'{group}vsHC']['p.value'].iloc[elec]
-    z_score = lmer_results[f'{group}vsHC']['z.ratio'].iloc[elec]
+    p_value = lmer_results['p.value'].iloc[elec]
+    if 'z.ratio' in lmer_results.columns:
+        statistic = lmer_results['z.ratio'].iloc[elec]
+    else:
+        statistic = lmer_results['t.ratio'].iloc[elec]
+    #statistic = lmer_results['z.ratio'].iloc[elec]
 
     if p_value < p_value_th:
-        data_lmer.append(z_score)
+        data_lmer.append(statistic)
     else:
         data_lmer.append(0)
 
@@ -51,7 +53,7 @@ def append_lmer_results(lmer_results, group, elec, p_value_th, data_lmer):
 if __name__ == "__main__":
     # Some parameters for the figure
     ncols = 6 
-    nrows = 5
+    nrows = 4  #5
 
     left = 0.06
     right = 0.11
@@ -67,7 +69,7 @@ def append_lmer_results(lmer_results, group, elec, p_value_th, data_lmer):
     spacing_y = 0.064 
 
     # Create figure
-    fig1 = plt.figure(figsize=(7.5, 5.5), dpi=300)
+    fig1 = plt.figure(figsize=(7.5, 5.5), dpi=150)
     current_bottom = bottom
 
     for row in range(2):
@@ -82,14 +84,18 @@ def append_lmer_results(lmer_results, group, elec, p_value_th, data_lmer):
         if row == 2 or row == 3:
             method = 'power_spectrum_parameterization_1'
         try:
-            data = pd.read_pickle(os.path.join('../data', method, 'emp_data_reduced.pkl'))
+            data = pd.read_pickle(os.path.join(results_path, method, 'emp_data_reduced.pkl'))
 
         except Exception as e:
             print(f'Error loading data for {method}: {e}')
             continue
 
+        data['Sensor'] = data['Sensor'].apply(lambda x: str(x))  # convert from np.str_ to str (to avoid warning when importing to R)
+        data = data[data['Group'] != 'MCI']  # remove MCI from dataset so p-value correction does not take it into account
+
         current_left = left
         for col in range(ncols):
+            print(f'\n----- Processing row {row}, column {col} -----\n')
             if col == 0 or col == 3:
                 group = 'ADMIL'
                 group_label = 'ADMIL'
@@ -125,37 +131,61 @@ def append_lmer_results(lmer_results, group, elec, p_value_th, data_lmer):
 
             if col >= 3:
                 var = 3 if row == 0 or row == 2 else 2
-
+
+            # Create a copy of the dataframe
+            data_preds = data.copy()
+            data_preds['Predictions'] = data_preds['Predictions'].apply(lambda x: x[var])
+
             # Statistical analysis
-            analysis = ncpi.Analysis(data)
+            analysis = ncpi.Analysis(data_preds)
             if statistical_analysis == 'lmer':
-                stat_results = analysis.lmer(control_group='HC', data_col='Predictions', data_index=var,
-                                        other_col=['ID', 'Group', 'Epoch', 'Sensor'],
-                                        models={
-                                            'mod00': 'Y ~ Group * Sensor + (1 | ID)',
-                                            'mod01': 'Y ~ Group * Sensor',
-                                            'mod02': 'Y ~ Group + Sensor + (1 | ID)',
-                                            'mod03': 'Y ~ Group + Sensor',
-                                        },
-                                        bic_models=["mod00", "mod01"],
-                                        anova_tests={
-                                            'test1': ["mod00", "mod01"],
-                                            'test2': ["mod02", "mod03"],
-                                        },
-                                        specs='~Group | Sensor')
+
+                # Reduce random effect structure with BIC and fixed effect structure with LRT (anova)
+                opt_f = analysis.lmer_selection(full_model='Predictions ~ Group * Sensor + (1 | ID)',
+                                                  numeric=['Predictions'],
+                                                  crit=None,
+                                                  fixed_crit='LRT',
+                                                  random_crit='BIC',
+                                                  include=['Sensor', 'Group'])
+                # Run post-hoc analyses using selected model
+                stat_results = analysis.lmer_tests(models=opt_f,
+                                              numeric=['Predictions'],
+                                              group_col='Group',
+                                              control_group='HC',
+                                              specs=['Group|Sensor'])
+
+                # opt_f_1 = 'Predictions ~ Group * Sensor + (1 | ID)'
+                # opt_f_2 = 'Predictions ~ Group + Sensor + (1 | ID)'
+                # stat_results = analysis.lmer_tests(models=[opt_f_1, opt_f_2],
+                #                               numeric=['Predictions'],
+                #                               group_col='Group',
+                #                               control_group='HC',
+                #                               specs=['Group|Sensor'])
+
+                # opt_f = 'Predictions ~ Group * Sensor + (1 | ID)'
+                # stat_results = analysis.lmer_tests(models=opt_f,
+                #                               numeric=['Predictions'],
+                #                               group_col='Group',
+                #                               control_group='HC',
+                #                               specs=['Group|Sensor'])
+
             elif statistical_analysis == 'cohend':
                 stat_results = analysis.cohend(control_group='HC', data_col='Predictions', data_index=var)
 
             data_stat = []
+
+            # Select the results for the current group
+            stat_results = stat_results['Group|Sensor'].query("contrast == @group + ' - HC'")
+
             # Extract sensor names
             empirical_sensors = data['Sensor'].unique()
             for elec in range(19):
                 # Find position of the electrode in the stat results
-                pos_results = np.where(stat_results[f'{group}vsHC']['Sensor'] == empirical_sensors[elec])[0]
+                pos_results = np.where(stat_results['Sensor'] == empirical_sensors[elec])[0]
 
                 if len(pos_results) > 0:
                     if statistical_analysis == 'lmer':
-                        data_stat = append_lmer_results(stat_results, group, pos_results[0], p_value_th, data_stat)
+                        data_stat = append_lmer_results(stat_results, pos_results[0], p_value_th, data_stat)
                     elif statistical_analysis == 'cohend':
                         data_stat.append(stat_results[f'{group}vsHC']['d'][pos_results[0]])
                 else:
@@ -238,6 +268,8 @@ def append_lmer_results(lmer_results, group, elec, p_value_th, data_lmer):
     fig1.add_artist(tauexc_line_f)
     fig1.add_artist(tauinh_line_f)
 
-    fig1.savefig('EEG_predictions.png')
+    # fig1.savefig('EEG_predictions.png')
+
+plt.show()
 
-    
+#plt.gcf().savefig(os.path.join(results_path, 'fig7.pdf'), bbox_inches='tight')

examples/statistics/lmer_example.py

@@ -0,0 +1,114 @@
+
+import numpy as np
+import pandas as pd
+
+from ncpi import Analysis
+
+####################################################################
+# Generate toy data
+
+np.random.seed(0)
+
+nsubj, nsensors, nepochs = 10, 3, 5
+id = np.repeat(np.arange(nsubj), nsensors*nepochs)
+sensor = np.tile(np.repeat(np.arange(nsensors), nepochs), nsubj)
+epoch = np.tile(np.arange(nepochs), nsensors*nsubj)
+data = pd.DataFrame([id, epoch, sensor], index=['id', 'epoch', 'sensor']).T
+group_size = nsubj//3
+data['gr'] = 'HC'
+data.loc[data['id'].isin(range(group_size)), 'gr'] = 'g1'
+data.loc[data['id'].isin(range(group_size, group_size * 2)), 'gr'] = 'g2'
+data['Y'] = np.random.normal(size=nsubj*nsensors*nepochs)
+# Group*Sensor effects
+halfsensor = nsensors//2
+data.loc[(data['gr'] == 'g1') & (data['sensor'] <= halfsensor), 'Y'] = data.loc[
+                    (data['gr'] == 'g1') & (data['sensor'] <= halfsensor), 'Y'] + np.random.normal(
+                    loc=10, scale=2, size=len((data.loc[(data['gr'] == 'g1') & (data['sensor'] <= halfsensor), 'Y'])))
+data.loc[data['gr'] == 'g2', 'Y'] = data.loc[data['gr'] == 'g2', 'Y'] + np.random.normal(
+                    loc=7, scale=.5, size = len(data.loc[data['gr'] == 'g2', 'Y']))
+# Make Y vary with epoch for HC and g2 but not for g1
+data.loc[data['gr'] != 'g1', 'Y'] = data.loc[data['gr'] != 'g1', 'Y'] + np.random.normal(
+                                loc=data.loc[data['gr'] != 'g1', 'epoch'], scale=.5)
+# Include subject-level deviations
+raneff = pd.DataFrame(np.random.normal(loc=0, scale=3, size=nsubj), index=range(nsubj), columns=['raneff'])
+data = data.join(raneff, on='id')
+data['Y'] = data['Y'] + np.random.normal(loc=data['raneff'], scale=.1)
+data.drop(columns='raneff', inplace=True)
+# Include covariates
+data['sex'] = np.random.choice(['F', 'M'], size=len(data))
+
+print('\nData:')
+print(data.head())
+
+####################################################################
+# Run analyses
+
+# Initialise class
+analysis = Analysis(data)
+
+# 1.
+print(f'\n\n{"Example 1":=^30}')
+# Reduce random effect structure with BIC
+print(f'\n{"Model selection":-^30}\n')
+opt_f = analysis.lmer_selection(full_model='Y ~ gr * epoch * sensor + sex + (1|id)',
+                                numeric=['epoch'],
+                                random_crit='BIC', fixed_crit=None)
+# Run post-hoc analyses using selected model
+print(f'\n{"Testing (1)":-^30}\n')
+results = analysis.lmer_tests(models=opt_f, group_col='gr', control_group='HC',
+                              numeric=['epoch'], specs=['epoch:gr', 'gr|sensor', 'x'])
+# Note: if a variable 'x' which is not in the models is given to specs, it is simply skipped.
+print(f'\n{"Testing (2)":-^30}\n')
+results = analysis.lmer_tests(models=opt_f, numeric=['epoch'], specs=['gr|sensor'])
+# Note: here we're not specifying the control group, so group comparisons include *all* combinations!
+
+# 2.
+print(f'\n\n{"Example 2":=^30}')
+# Reduce random effect structure with BIC and fixed effect structure with LRT (anova)
+print(f'\n{"Model selection":-^30}\n')
+opt_f = analysis.lmer_selection(full_model='Y ~ gr * epoch * sensor + sex + (1|id)',
+                                numeric=['epoch'],
+                                random_crit='BIC', fixed_crit='LRT')
+# Run post-hoc analyses using selected model
+print(f'\n{"Testing":-^30}\n')
+results = analysis.lmer_tests(models=opt_f, group_col='gr', control_group='HC',
+                              numeric=['epoch'], specs=['gr|sensor', 'epoch:gr', 'sex'])
+# Note: as sex is not present in the model anymore, the test is automatically skipped.
+
+# 3.
+print(f'\n\n{"Example 3":=^30}')
+# Reduce random effect structure with BIC and fixed effect structure with LRT (anova),
+# but always keeping the sensor term in the model
+print(f'\n{"Model selection":-^30}\n')
+opt_f = analysis.lmer_selection(full_model='Y ~ gr * epoch * sensor + sex + (1|id)',
+                                numeric=['epoch'],
+                                random_crit='BIC', fixed_crit='LRT', include=['sensor'])
+# Run post-hoc analyses using selected model
+print(f'\n{"Testing":-^30}\n')
+results = analysis.lmer_tests(models=opt_f, group_col='gr',
+                              control_group='HC', numeric=['epoch'])
+# Note: no specs given, so group_col is used.
+
+# 4.
+print(f'\n\n{"Example 4":=^30}')
+# Use BIC to compare only a specified set of possible models
+# (i.e., without exploring all their subsets via backward selection),
+# then run post-hoc analyses
+results = analysis.lmer_tests(models=['Y ~ gr * epoch * sensor + sex + (1|sensor) + (1|id)',
+                                      'Y ~ gr * epoch * sensor + (1|sex) + (1|id)'],  # Best model is selected using BIC
+                group_col='gr', control_group='HC', numeric=['epoch'],
+                specs=['gr|sensor', 'epoch:gr'])
+
+# 5.
+print(f'\n\n{"Example 5":=^30}')
+# Reduce random effect structure with BIC and fixed effect structure with LRT (anova)
+print(f'\n{"Model selection":-^30}\n')
+opt_f = analysis.lmer_selection(full_model='Y ~ gr * epoch * sensor + sex + (1|id)',
+                                numeric=['epoch'],
+                                random_crit='BIC', fixed_crit='LRT')
+# Run post-hoc analyses using selected model
+print(f'\n{"Testing":-^30}\n')
+results = analysis.lmer_tests(models=opt_f, group_col='gr', control_group='HC',
+                              numeric=['epoch'], specs=['sex'])
+# Note: as sex is not present in the model, the test is automatically skipped.
+# No specs left: empty output and warning message.