<jupyter_start><jupyter_text>Update February 13th 2019--We performed all classification experiments again the evoked response, but we selected components using a Kbest with 25 electrodes before training the classifiers.  We train classifiers on one condition, and test it on the same or other condition (eg. Train Regular, Test Regular), so this is why we get four graphs. One interesting part would be to observe a group difference in the respect to which classifiers in one group are able to generalize across conditions, or not. For instance, one interesting hypothesis to test would be : - Classifiers trained on controls do generalize from the pattern trained in Regular when tested in Irregular, or the opposite. - Classifiers trained on Patients data are less/not able to generalize across conditions.<jupyter_code>import numpy as np
from matplotlib import pyplot as plt 
%matplotlib inline 
import os 
from statsmodels.stats import multitest

statspath = '/home/nfarrugi/datasets/mpi_pd_cueing/results_fkbest25/'<jupyter_output><empty_output><jupyter_text>Step 0 - preparing data  fields for the saved stats files are : - scores- proba- pval- UFor each of them, layout is 4 (Regreg,RegIrreg,IrregReg,IrregIrreg) x CV x TestTime x TrainTime Group files per subject type (patients or control)<jupyter_code>allsubj = os.listdir(statspath)

controls_id = []
patients_id = []

for cursubj in allsubj:
    if (cursubj[9]=='k'):
        controls_id.append(cursubj)
    else:
        if cursubj[10]=='1':
            patients_id.append(cursubj)
            
len(patients_id)
len(controls_id)
def diagonal_group(subjlist,statspath,thresh_unc = 0.01,timevec=None,showall=False,threshold_zeros=40,doplots=True,keepall=True,onlyregreg=False):


    if (showall):
        print("Showing all subjects")
    else:
        print("Showing subjects with less than %d FDR significant points" % threshold_zeros)
    
    diag_sig_fdr_all = []
    diag_sig_unc_all = []
    diag_ROC_all = []
    for subj in subjlist:
        cursubj = os.path.join(statspath,subj)
        curroc = np.load(cursubj)['scores'].mean(axis=1)
        curpvals = np.load(cursubj)['pval'].mean(axis=1) # Calculating mean p-value accross CV splits 
        diag_pvals = np.stack([np.diag(curmat) for curmat in curpvals]) # extract diagonal classification for each condition pair

        diag_ROC = np.stack([np.diag(curmat) for curmat in curroc]) # extract diagonal classification
        
        diag_sig_fdr = []
        diag_sig_unc = []
        
        if onlyregreg:
            H,p,_,_=multitest.multipletests(diag_pvals[2],method='fdr_bh')
            H_unc = diag_pvals[2] < thresh_unc
            diag_sig_unc.append(H_unc)
            diag_sig_fdr.append(H)        
        else:
            for curdiag in diag_pvals:
                H,p,_,_=multitest.multipletests(curdiag,method='fdr_bh')
                H_unc = curdiag < thresh_unc
                diag_sig_unc.append(H_unc)
                diag_sig_fdr.append(H)

        
        nonzeros = np.sum(np.stack(diag_sig_fdr))
        if showall:
            print("Subj %s number of non-zeros FDR : %d " %(subj,nonzeros))
                        
        elif nonzeros < threshold_zeros:
            print("Subj %s has zero or few significant FDR : %d " %(subj,nonzeros))
            
            if (nonzeros >= 0 ) and doplots:
                #plt.subplot(1,2,1)
                plt.plot(timevec,diag_sig_fdr[0])
                #plt.subplot(1,2,2)
                plt.plot(timevec,diag_ROC[0])
                plt.show()
            
        
        #print("Subj %s number of non-zeros Uncorrected : %d " %(subj,np.sum(np.stack(diag_sig_unc))))
        
        if (nonzeros > threshold_zeros) or keepall:
            diag_sig_fdr_all.append(np.stack(diag_sig_fdr))
            diag_sig_unc_all.append(np.stack(diag_sig_unc))
            diag_ROC_all.append(diag_ROC)

    
    return np.sum(np.stack(diag_sig_fdr_all),axis=0),np.sum(np.stack(diag_sig_unc_all),axis=0),np.stack(diag_ROC_all)
def plot_ROC_allcond(allscores,timepoints,figtitle='Default Title'):

    fig, (ax) = plt.subplots(ncols=2,nrows=2,figsize=(10,10))

    titles = ['Train Regular Test Regular',
              'Train Regular Test Irregular',
              'Train Irregular Test Regular',
              'Train Irregular Test Irregular']

    for i in range(4):
        curax = ax.ravel()[i]

        im = curax.matshow(allscores.mean(axis=1)[i], cmap='RdBu_r', vmin=0.2,vmax=0.8,origin='lower',
                            extent=timepoints)
        curax.axhline(0., color='k')
        curax.axvline(0., color='k')
        curax.xaxis.set_ticks_position('bottom')
        curax.set_xlabel('Testing Time (s)')
        curax.set_ylabel('Training Time (s)')
        curax.set_title(titles[i])
        plt.colorbar(im, ax=curax)
    fig.suptitle(figtitle)
    return fig
def plot_ROC_group(subjlist,statspath):

    diag_sig_fdr_all = []
    diag_sig_unc_all = []
    diag_ROC_all = []
    
    timepoints = np.array([-0.05      ,  0.52000002, -0.05      ,  0.52000002])
    for subj in subjlist:
        cursubj = os.path.join(statspath,subj)
        curroc = np.load(cursubj)['scores']
        curpvals = np.load(cursubj)['pval']
        #print(subj)
        
        plot_ROC_allcond(curroc,timepoints,subj)<jupyter_output><empty_output><jupyter_text>Part1-- Individual plots of ROC curves for temporal generalization using evoked resp + PCA80perc to train classifiers (averaged over four consecutive trials)--<jupyter_code>#plot_ROC_group(controls_id,statspath)
#plot_ROC_group(patients_id,statspath)<jupyter_output><empty_output><jupyter_text>Individual ROC curves for temporal generalization on evoked response to train classifiers without PCA ( just showing again previous results here to compare)--<jupyter_code>#plot_ROC_group(avcontrols_id,avstatspath)
#plot_ROC_group(avpatients_id,avstatspath)
def four_plots(timevec,subtitles,diag_sig,mytitle):
    plt.figure(figsize=(10,10))

    for i in range(1,5):
        plt.subplot(2,2,i)
        plt.plot(timevec,diag_sig[i-1])
        plt.title(subtitles[i-1])
    plt.suptitle('Diagonal decoding - %s' % mytitle)
    plt.show()
def four_errorplots(timevec,subtitles,diag_roc,mytitle):
    plt.figure(figsize=(10,10))

    for i in range(1,5):
        plt.subplot(2,2,i)
        plt.errorbar(timevec,diag_roc[:,i-1,:].mean(axis=0),yerr=diag_roc[:,i-1,:].std(axis=0),errorevery=5)
        plt.title(subtitles[i-1])
        plt.grid(True)
    plt.suptitle('Diagonal decoding - %s' % mytitle)
    plt.show()<jupyter_output><empty_output><jupyter_text>Part2-- Summarizing statistics using the beta results-- Step 1 - Analysis of "Diagonal" ClassificationWe average the diagonal of the p values of the within-subjects Mann Whittney U-Test, and average them across CV folds. We correct the resulting p values using FDR Benjamin Hochberg correction, and sum the resulting vectors across all subjects. We also compute an uncorrected version with a threshold of $p < 0.01$<jupyter_code>diag_sig_fdr_controls,diag_sig_unc_controls,diag_ROC_controls = diagonal_group(controls_id,statspath,0.05,doplots=False)
diag_sig_fdr_patients,diag_sig_unc_patients,diag_ROC_patients = diagonal_group(patients_id,statspath,0.05,doplots=False)<jupyter_output>Showing subjects with less than 40 FDR significant points
Subj 180802_07k1_formal_conditionwise.npz has zero or few significant FDR : 23 
Subj 180802_02k1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802_11k1_formal_conditionwise.npz has zero or few significant FDR : 18 
Subj 180802_15k1_formal_conditionwise.npz has zero or few significant FDR : 35 
Subj 180802_18k1_formal_conditionwise.npz has zero or few significant FDR : 26 
Showing subjects with less than 40 FDR significant points
Subj 180802_12c1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802_15c1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802_17c1_formal_conditionwise.npz has zero or few significant FDR : 4 
Subj 180802_11c1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802_14c1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802_04c1_formal_conditionwise.npz has zero or few significant FDR : 1 
Subj 180802[...]<jupyter_text>5 controls and 11 patients don't have sig FDR, which is a bit better than when taking all electrodes.  Uncomment the following cell to output number of FDR corrected <NAME> whithin subjects statistics for all subjects (to give an idea how many signficant values we get for the ones who work well)<jupyter_code>diag_sig_fdr_controls,diag_sig_unc_controls,_ = diagonal_group(controls_id,statspath,0.05,showall=True)
diag_sig_fdr_patients,diag_sig_unc_patients,_ = diagonal_group(patients_id,statspath,0.05,showall=True)<jupyter_output>Showing all subjects
Subj 180802_07k1_formal_conditionwise.npz number of non-zeros FDR : 23 
Subj 180802_04k1_formal_conditionwise.npz number of non-zeros FDR : 227 
Subj 180802_13k1_formal_conditionwise.npz number of non-zeros FDR : 357 
Subj 180802_05k1_formal_conditionwise.npz number of non-zeros FDR : 140 
Subj 180802_06k1_formal_conditionwise.npz number of non-zeros FDR : 53 
Subj 180802_20k1_formal_conditionwise.npz number of non-zeros FDR : 286 
Subj 180802_09k1_formal_conditionwise.npz number of non-zeros FDR : 119 
Subj 180802_08k1_formal_conditionwise.npz number of non-zeros FDR : 245 
Subj 180802_02k1_formal_conditionwise.npz number of non-zeros FDR : 0 
Subj 180802_12k1_formal_conditionwise.npz number of non-zeros FDR : 143 
Subj 180802_11k1_formal_conditionwise.npz number of non-zeros FDR : 18 
Subj 180802_01k1_formal_conditionwise.npz number of non-zeros FDR : 387 
Subj 180802_14k1_formal_conditionwise.npz number of non-zeros FDR : 91 
Subj 180802_15k1_formal_conditionwis[...]<jupyter_text>bottomline : Kbest DID boost SNR ! Maybe 25 electrodes is too much / too less ? Let's try with 40 electrodes.<jupyter_code>tmin = -0.05
tmax = 0.52

timevec = np.linspace(tmin,tmax,diag_sig_fdr_controls.shape[1])
subtitles = ['Train Regular Test Regular',
            'Train Regular Test Irregular',
            'Train Irregular Test Regular',
            'Train Irregular Test Irregular']<jupyter_output><empty_output><jupyter_text>Uncomment the following cell to just plot the sum of significant Mann Wittney subjects as a function of time (FDR corrected)<jupyter_code>##four_plots(timevec,subtitles,diag_sig_fdr_controls,'Controls - FDR')
##four_plots(timevec,subtitles,diag_sig_fdr_patients,'Patients')<jupyter_output><empty_output><jupyter_text>Uncomment the following lines to look at the same counts without FDR correction<jupyter_code>four_plots(timevec,subtitles,diag_sig_unc_controls,'Controls-unc')
four_plots(timevec,subtitles,diag_sig_unc_patients,'Patients unc')<jupyter_output><empty_output><jupyter_text>We can examine the ROC AUC scores for the same subjects, to check whether this phenomenon has anything to do with classification performance (it should!). For that, we overlay the FDR-corrected timewise significance of the Mann-Whitney U-test with the ROC AUC performance of the problematic subjects<jupyter_code>#diag_sig_fdr_controls,diag_sig_unc_controls,diag_ROC_controls = diagonal_group(controls_id,statspath,0.05,timevec=timevec,doplots=True)
#diag_sig_fdr_patients,diag_sig_unc_patients,diag_ROC_patients = diagonal_group(patients_id,statspath,0.05,timevec=timevec,doplots=True)<jupyter_output><empty_output><jupyter_text>To illustrate this at the group level, let's now plot the average and std of ROC AUC for the whole groups.<jupyter_code>four_errorplots(timevec,subtitles,diag_roc=diag_ROC_controls,mytitle='Controls AUC')<jupyter_output><empty_output><jupyter_text>For patients :<jupyter_code>four_errorplots(timevec,subtitles,diag_roc=diag_ROC_patients,mytitle='Patients - AUC')<jupyter_output><empty_output><jupyter_text>Let's have a look at the averages only, but superimposing patients and controls on the same plots, zooming in on the region 350ms to 500 ms. We will superimpose the result of an independent t-test at each time point, between groups<jupyter_code>from scipy.stats import ttest_ind

def twogroups_fourplots(timevec,subtitles,diag_roc1,diag_roc2,mytitle,tstart=None,tend=None,ttest=False):
    plt.figure(figsize=(10,10))

    ngroup1 = diag_roc1.shape[0]
    ngroup2 = diag_roc2.shape[0]
    
    if tstart is None:
        tstart=0
    if tend is None:
        tend=timevec.shape[0]
    
    for i in range(1,5):
        plt.subplot(2,2,i)
        plt.plot(timevec[tstart:tend],diag_roc1[:,i-1,tstart:tend].mean(axis=0))
        plt.plot(timevec[tstart:tend],diag_roc2[:,i-1,tstart:tend].mean(axis=0))
        
        if ttest:
            
            
            tstat,p_unc = ttest_ind(a=diag_roc1[:,i-1,tstart:tend],b=diag_roc2[:,i-1,tstart:tend],equal_var=True)
            
            timevec_sub = timevec[tstart:tend]
            H = (p_unc<0.05)
            sig_ind = timevec_sub[H]
            sig= np.ones_like(sig_ind)*0.9
            
            #FDR correction            
            #H,p_fdr,_,_=multitest.multipletests(p_unc,method='fdr_bh')
            
            plt.plot(sig_ind,sig,'*')
            
            
        plt.title(subtitles[i-1])
        plt.grid(True)
    plt.legend(['controls (n=%d)' % ngroup1,'patients (n=%d)' % ngroup2,'ttest uncorr.'],loc='center left', bbox_to_anchor=(1, 0.5))
    
    plt.suptitle('Diagonal decoding - %s' % mytitle)
    plt.show()
twogroups_fourplots(timevec,subtitles,diag_ROC_controls,diag_ROC_patients,'Comparison - All subjects',tstart=120,tend=250,ttest=True)<jupyter_output><empty_output><jupyter_text>Part3-- Evoked responses-- High SNR subjects<jupyter_code>diag_sig_fdr_controls,diag_sig_unc_controls,diag_ROC_controls = diagonal_group(controls_id,statspath,0.05,doplots=False,keepall=False)
diag_sig_fdr_patients,diag_sig_unc_patients,diag_ROC_patients = diagonal_group(patients_id,statspath,0.05,doplots=False,keepall=False)

four_errorplots(timevec,subtitles,diag_roc=diag_ROC_controls,mytitle='Controls AUC')
four_errorplots(timevec,subtitles,diag_roc=diag_ROC_patients,mytitle='Patients - AUC')

twogroups_fourplots(timevec,subtitles,diag_ROC_controls,diag_ROC_patients,'Comparison - High SnR subjects',tstart=120,tend=250,ttest=True)
diag_sig_fdr_controls,diag_sig_unc_controls,diag_ROC_controls = diagonal_group(controls_id,statspath,0.05,doplots=False,keepall=True)
diag_sig_fdr_patients,diag_sig_unc_patients,diag_ROC_patients = diagonal_group(patients_id,statspath,0.05,doplots=False,keepall=True)

four_errorplots(timevec,subtitles,diag_roc=diag_ROC_controls,mytitle='Controls AUC')
four_errorplots(timevec,subtitles,diag_roc=diag_ROC_patients,mytitle='Patients - AUC')

twogroups_fourplots(timevec,subtitles,diag_ROC_controls,diag_ROC_patients,'Comparison - All subjects',tstart=120,tend=250,ttest=True)<jupyter_output>Showing subjects with less than 40 FDR significant points
Subj 180802_07k1_formal_conditionwise.npz has zero or few significant FDR : 18 
Subj 180802_04k1_formal_conditionwise.npz has zero or few significant FDR : 35 
Subj 180802_06k1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802_02k1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802_12k1_formal_conditionwise.npz has zero or few significant FDR : 21 
Subj 180802_11k1_formal_conditionwise.npz has zero or few significant FDR : 6 
Subj 180802_15k1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802_17k1_formal_conditionwise.npz has zero or few significant FDR : 39 
Subj 180802_16k1_formal_conditionwise.npz has zero or few significant FDR : 6 
Subj 180802_18k1_formal_conditionwise.npz has zero or few significant FDR : 0 
Showing subjects with less than 40 FDR significant points
Subj 180802_12c1_formal_conditionwise.npz has zero or few significant FDR : 0 
Subj 180802[...]
{'fraction_non_alphanumeric': 0.08879824286162535, 'fraction_numerical': 0.03897081895199247, 'mean_word_length': 4.409368635437882, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '242065', 'license': 'Apache-2.0', 'n_tokens_mistral': 5992, 'n_tokens_neox': 5361, 'n_words': 1198, 'path': 'MNE-TF-Statistics-Averages-Kbest25.ipynb', 'repo_name': 'nicofarr/pd_osc_decoding', 'size': '15941'}