notebook.community

Edit and run 





In [1]:



    


%matplotlib inline
main_path = '../../../Dropbox/McGill-publication/Papers/Multisite/'
path_out = main_path + 'figures/'



    











In [2]:



    


from matplotlib.pylab import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.io
from mpl_toolkits.axes_grid1 import ImageGrid
import matplotlib.gridspec as gridspec

def errorfill(x, y, yerr, color=None, alpha_fill=0.3, ax=None, label=''):
    ax = ax if ax is not None else plt.gca()
    if color is None:
        color = ax._get_lines.color_cycle.next()
    if np.isscalar(yerr) or len(yerr) == len(y):
        ymin = y - yerr
        ymax = y + yerr
    elif len(yerr) == 2:
        ymin, ymax = yerr
    ax.plot(x, y, color=color, label=label, lw=2)
    ax.fill_between(x, ymax, ymin, color=color, alpha=alpha_fill)
    
def load_file(path):
    mat = scipy.io.loadmat(path)
    return mat

def gen_fig(path_in, file_in,path_out, file_out, conn = 0, pvalidx=0):
    labels_box = ['H0 multi','H0 Dummy', 'H0 METAL']
    lw = 2
    mat = load_file(path_in + file_in)
    pvals = mat.get('list_p')[0]
    fig = plt.figure(figsize=(8, 4))
    sns.set_context('paper',font_scale=1.7)
    plt.subplot(1, 2, 1)
    #ax.set_xticks(np.arange(0,1,0.1))
    #ax.set_yticks(np.arange(0,1.,0.1))
    plt.xlim(0,1.5)
    effect_size = np.arange(0.0, 1.501, 0.01)
    plt.ylabel('Detection power')
    plt.xlabel('Effect size (Cohen\'s ' + r'$d$)')
    plt.plot(effect_size, mat.get('sens_monosite')[pvalidx,:,conn],label="Monosite", lw=lw, c='r')
    #plt.plot(effect_size, mat.get('sens_multisite_nocorr')[pvalidx,:,conn],label="Multisite no correction",lw=lw, c='k')
    plt.plot(effect_size, mat.get('sens_multisite_dummyvar')[pvalidx,:,conn],label="Multisite (7 sites)",lw=lw, c='b')
    #plt.plot(effect_size, mat.get('sens_multisite_metal')[pvalidx,:,conn],label="Multisite METAL",lw=lw, c='m')
    plt.grid(True)
    plt.title('p-value ' +  str(pvals[pvalidx]))
    
    # fully use the given bounding box.
    #plt.legend(bbox_to_anchor=(1.05, 1), loc=4, mode="expand", borderaxespad=0.)
    plt.legend( bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0., frameon=False)
    # do the box plot
    plt.subplot(2, 2, 4)
    ax = fig.gca()
    #xtickNames = plt.setp(ax, xticklabels=np.repeat(labels_box, 1))
    #plt.setp(xtickNames, rotation=45, fontsize=12)
    plt.xlabel(r'Effect size (Cohen\'s $d$)')
    #h0_metal =  mat.get('sens_multisite_h0_metal')[pvalidx,:,conn]
    #h0_multi  = mat.get('sens_multisite_h0')[pvalidx,:,conn]
    h0_dummy  = np.mean(mat.get('sens_multisite_h0_dummy')[pvalidx,:,conn],axis=0)
    #boxplot((h0_multi,h0_dummy,h0_metal))
    boxplot((h0_dummy))
    ax.set_xticklabels(labels_box, rotation=0, ha='center')
    ax.set_yscale('log')
    ax.set_ylim((10**(-3),10**(0)))
    plt.ylabel('Detection power (log)')
    plt.title('H0 hypotheses')
    plt.tight_layout()
    #plt.savefig(path_out + file_out)
    #plt.show()

def gen_fig_allconn(path_in, file_in, path_out, file_out, conn=0, pvalidx=0, efflim=1.5,nsites=7):
    labels_box = ['H0 multi','H0 Dummy', 'H0 METAL']
    lw = 2
    mat = load_file(path_in + file_in)
    pvals = mat.get('list_p')[0]
    
    fig = plt.figure(figsize=(5,5))
    sns.set_context('paper',font_scale=1.7)
    plt.subplot2grid((5,2), (1,0),colspan=2,rowspan=4)
    ax = fig.gca()
    #ax.set_xticks(np.arange(0,1,0.1))
    #ax.set_yticks(np.arange(0,1.,0.1))
    if efflim==2:
        plt.xlim(0,2)
        effect_size = np.arange(0.0, 2.001, 0.01)
    else:
        plt.xlim(0,1.5)
        effect_size = np.arange(0.0, 1.501, 0.01)
    
    plt.ylabel('Detection power')
    plt.xlabel('Effect size (Cohen\'s ' + r'$d$)')
    ax.set_ylim((0,1))
    tmpdata = mat.get('sens_monosite')[pvalidx,:,:]
    errorfill(effect_size, tmpdata.mean(axis=1)[:len(effect_size)],tmpdata.std(axis=1)[:len(effect_size)],label="Monosite", color='r')
    #tmpdata = mat.get('sens_multisite_nocorr')[pvalidx,:,:]
    #errorfill(effect_size, tmpdata.mean(axis=1),tmpdata.std(axis=1),label="Multisite no correction", color='k')
    tmpdata = mat.get('sens_multisite_dummyvar')[pvalidx,:,:]
    errorfill(effect_size, tmpdata.mean(axis=1)[:len(effect_size)],tmpdata.std(axis=1)[:len(effect_size)],label="Multisite ("+str(nsites)+" sites)", color='b')
    #tmpdata = mat.get('sens_multisite_metal')[pvalidx,:,:]
    #errorfill(effect_size, tmpdata.mean(axis=1),tmpdata.std(axis=1),label="Multisite METAL", color='m')
    
    plt.grid(True)
    plt.title(r'False-positive rate $\alpha=$' +  str(pvals[pvalidx]))
    
    # fully use the given bounding box.
    plt.legend(bbox_to_anchor=(0., 1.10, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.,frameon=False)
    #plt.legend( bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0., frameon=False)
    # do the box plot
    #plt.subplot2grid((7,2), (5,0),colspan=2,rowspan=2)
    #ax = fig.gca()
    ##xtickNames = plt.setp(ax, xticklabels=np.repeat(labels_box, 1))
    ##plt.setp(xtickNames, rotation=45, fontsize=12)
    #plt.xlabel('Effect size (cohen\'s d)')
    #h0_metal =  np.ndarray.flatten(np.mean(mat.get('sens_multisite_h0_metal')[pvalidx,:,:],axis=0))
    #h0_multi  = np.ndarray.flatten(np.mean(mat.get('sens_multisite_h0')[pvalidx,:,:],axis=0))
    #h0_dummy  = np.ndarray.flatten(np.mean(mat.get('sens_multisite_h0_dummy')[pvalidx,:,:],axis=0))
    ##boxplot((h0_multi,h0_dummy,h0_metal))
    #boxplot((h0_dummy))
    #ax.set_xticklabels(labels_box, rotation=0, ha='center')
    #ax.set_yscale('log')
    #ax.set_ylim((10**(-5),10**(0)))
    #plt.ylabel('Detection power (log)')
    #plt.title('H0 hypotheses')
    #plt.grid(True)
    print('monosite h0: ',np.mean(mat.get('sens_monosite_h0')[pvalidx,:,:]))
    plt.tight_layout()
    plt.savefig(path_out + file_out,dpi=150)
    

def gen_fig_full_sim(path_in, file_in,path_out, file_out):
    labels_box = ['H0 multi','H0 Dummy', 'H0 METAL']
    lw = 2
    mat = load_file(path_in + file_in)
    pvals = 0.001
    #fig = plt.figure(figsize=(5.5,7))
    fig = plt.figure(figsize=(5,5))
    sns.set_context('paper',font_scale=1.7)
    plt.subplot2grid((5,2), (1,0),colspan=2,rowspan=4)
    ax = fig.gca()
    plt.xlim(0,1.5)
    effect_size = np.arange(0.0, 1.501, 0.01)
    ax.set_ylim((0,1))
    plt.ylabel('Detection power')
    plt.xlabel('Effect size (Cohen\'s ' + r'$d$)')
    plt.plot(effect_size, mat.get('sens_mono'),label="Monosite",lw=lw, color='r')
    #plt.plot(effect_size, mat.get('sens_multi'),label="Multisite no correction",lw=lw, c='k')
    plt.plot(effect_size, mat.get('sens_multi_dummy'),label="Multisite (7 sites)",lw=lw, c='b')
    #plt.plot(effect_size, mat.get('sens_multi_metal'),label="Multisite METAL",lw=lw, c='m')
    plt.grid(True)
    plt.title(r'False-positive rate $\alpha=$' +  str(pvals))
    
    # fully use the given bounding box.
    plt.legend(bbox_to_anchor=(0., 1.10, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.,frameon=False)
    ##plt.legend( bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0., frameon=False)
    ## do the box plot
    #plt.subplot2grid((7,2), (5,0),colspan=2,rowspan=2)
    #ax = plt.gca()
    ##ax = fig.gca()
    ##xtickNames = plt.setp(ax, xticklabels=np.repeat(labels_box, 1))
    ##plt.setp(xtickNames, rotation=45, fontsize=12)
    #plt.xlabel('Effect size (cohen\'s d)')
    ##h0_metal =  np.mean(mat.get('sens_multi_h0_metal'),axis=0)
    ##h0_multi  = np.mean(mat.get('sens_multi_h0'),axis=0)
    ##h0_dummy  = np.mean(mat.get('sens_multi_h0_dummy'),axis=0)
    ##boxplot((h0_multi,h0_dummy,h0_metal))
    #boxplot((h0_dummy))
    #ax.set_xticklabels(labels_box, rotation=0, ha='center')
    #ax.set_yscale('log')
    #ax.set_ylim((10**(-5),10**(0)))
    #plt.ylabel('Detection power (log)')
    #plt.title('H0 hypotheses')
    #plt.grid(True)
    print('monosite h0: ',np.mean(mat.get('sens_monosite_h0'),axis=0))
    plt.tight_layout()
    plt.savefig(path_out + file_out,dpi=150)
    ##plt.show()
    
def geteffectsize(path_in, file_in,mono=True,pidx=0):
    point_list = []
    mat = load_file(path_in + file_in)
    pvals = mat.get('list_p')[0]
    #print pvals
    if mono:
        pow_mask = mat.get('sens_monosite')[pidx]>=0.8 
    else:
        pow_mask = mat.get('sens_multisite_dummyvar')[pidx]>=0.8 
    
    for i in range(pow_mask.shape[1]):
        #print np.where(pow_mask[:,i])[0][0]
        point_list.append(mat.get('list_effect_size')[0][np.where(pow_mask[:,i])[0][0]])
    
    return np.mean(point_list),np.std(point_list)
#sens_monosite 
#sens_multisite_dummyvar



    











In [3]:



    


from scipy.interpolate import interp1d

import statsmodels.stats.power as pow
import pandas as pan
import math as math

plt.figure()
x = np.linspace(0, 2 * np.pi)
y_sin = np.sin(x)
y_cos = np.cos(x)
print(x.shape,y_cos.shape)
errorfill(x, y_sin, 0.5,'y')
errorfill(x, y_cos, 0.2,'c')

a = np.array([[1,2], [3,4]])
print a
print np.ndarray.flatten(a)



    


















    






((50,), (50,))
[[1 2]
 [3 4]]
[1 2 3 4]

alpha at 0.001 80% detection power





In [5]:



    


import seaborn; seaborn.set()
import seaborn as sns
sns.set_style("whitegrid")
sns.set_color_codes("muted")
 

path_tmp = main_path + 'figures/default_sim/samps_s40_50pct/'
m1,s1 = geteffectsize(path_tmp,'results_simu_power.mat',True)
path_tmp = main_path + 'figures/default_sim/samps_s80_50pct/'
m2,s2 = geteffectsize(path_tmp,'results_simu_power.mat',True) 
path_tmp = main_path + 'figures/default_sim/samps_s120_50pct/'
m3,s3 = geteffectsize(path_tmp,'results_simu_power.mat',True) 

path_tmp = main_path + 'figures/default_sim/samps_s40_50pct/'
mm1,sm1 = geteffectsize(path_tmp,'results_simu_power.mat',False) 
path_tmp = main_path + 'figures/default_sim/samps_s80_50pct/'
mm2,sm2 = geteffectsize(path_tmp,'results_simu_power.mat',False) 
path_tmp = main_path + 'figures/default_sim/samps_s120_50pct/'
mm3,sm3 = geteffectsize(path_tmp,'results_simu_power.mat',False)


yref=[]
x = np.linspace(40, 120, 3)
y = np.array([m1,m2,m3])
y2 = np.array([mm1,mm2,mm3])

f1 = interp1d(x, y, kind='linear')
f2 = interp1d(x, y2, kind='linear')

xnew = np.linspace(40, 120, 40)

for i in xnew:
    yref.append((pow.TTestIndPower().solve_power(nobs1=i/2, power=0.8, ratio=1, alpha=0.001)))
yref  = np.array(yref)

plt.figure(figsize=(6,6))
sns.set_context('paper',font_scale=1.8)
plt.plot(xnew,yref,'k--',lw=2)

#plt.scatter(x,y,c="r",s=50)
#plt.scatter(x[1:],y2,c="b",s=50)
plt.errorbar(x, y, yerr=[s1,s2,s3], fmt='rD')
plt.errorbar(x, y2, yerr=[sm1,sm2,sm3], fmt='bD')

plt.legend(['parametric','monosite','multisite (7 sites)'],bbox_to_anchor=(0., 1.05, 1., .002),handletextpad=0.01, loc=3, ncol=3, mode="expand",borderaxespad=0.,frameon=False)
plt.xlabel('Total # subjects')
plt.ylabel('Effect size (Cohen\'s ' + r'$d$)')
plt.grid(True)
axes = plt.gca()
axes.set_xlim([20,140])
axes.set_ylim([0,1.8])

plt.title(r'False-positive rate $\alpha=0.001$')
plt.savefig(path_out + 'samplesize_effectsize_pow80_alpha001.png',dpi=150)



    


















    
























In [4]:



    


import seaborn; seaborn.set()
import seaborn as sns
sns.set_style("whitegrid")
sns.set_color_codes("muted")
 
mat = load_file(main_path + 'figures/default_sim/samps_s120_15pct/results_simu_power.mat')



    


















    






/home/cdansereau/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))

















In [86]:



    


ls ../../../Dropbox/McGill-publication/Papers/Multisite/figures/default_sim/



    


















    






logs/                         samps_s33_25_75/  samps_s40_50pct/
n_subject_estimation_bis.mat  samps_s33_30pct/  samps_s40_rnd/
samps_s120_15pct/             samps_s33_50pct/  samps_s80_15pct/
samps_s120_30pct/             samps_s33_75_25/  samps_s80_30pct/
samps_s120_50pct/             samps_s33_rnd/    samps_s80_50pct/
samps_s120_rnd/               samps_s40_15pct/  samps_s80_rnd/
samps_s33_15pct/              samps_s40_30pct/

















In [81]:



    


print mat['sens_multisite_h0_dummy'][2,:,:]

mat.keys()



    


















    






[[ 0.053  0.041  0.053 ...,  0.062  0.051  0.053]
 [ 0.046  0.054  0.059 ...,  0.053  0.045  0.045]
 [ 0.051  0.05   0.042 ...,  0.042  0.055  0.052]
 ..., 
 [ 0.054  0.055  0.048 ...,  0.057  0.042  0.055]
 [ 0.06   0.039  0.045 ...,  0.056  0.061  0.047]
 [ 0.05   0.044  0.045 ...,  0.059  0.04   0.052]]










    
Out[81]:








['list_effect_size',
 'sens_monosite_h0',
 'sens_multisite_dummyvar',
 'sens_multisite_h0_dummy',
 'sens_multisite_h0',
 'sens_multisite_nocorr',
 '__header__',
 'sens_multisite_metal',
 '__globals__',
 'sens_monosite',
 'sens_multisite_h0_metal',
 'list_p',
 '__version__']

















In [90]:



    


def plot_boxplot_h0(mat,path_out,N,W,alpha_idx=0):
    alpha = [0.001,0.01,0.05]
    #fig, ax  = plt.figure(figsize=(4,4))
    fig, ax = plt.subplots(figsize=(4, 4))
    plt.plot(alpha[alpha_idx]*np.ones((8,1)),'k--')
    width_=1
    data = [mat['sens_monosite_h0'][alpha_idx,:,:].mean(),mat['sens_multisite_h0_dummy'][alpha_idx,:,:].mean(),mat['sens_multisite_h0'][alpha_idx,:,:].mean()]
    #plt.boxplot([mat['sens_monosite_h0'][alpha_idx,:,:],mat['sens_multisite_h0_dummy'][alpha_idx,:,:],mat['sens_multisite_h0'][alpha_idx,:,:]],notch=1,labels=('Monosite','Multisite\ncorrection','Multisite\nno correction'))
    plt.bar([1,3,5],data,width_)
    ax.set_xticks((1.5,3.5,5.5))
    ax.set_xticklabels(('Monosite','Multisite\ncorrection','Multisite\nno correction'))
    ax.set_ylabel('Detection power')
    ax.set_xlabel('N='+str(N)+', W='+str(W)+'%')
    n = mat['sens_monosite_h0'][alpha_idx,:,:].reshape(-1,1).shape[0]
    ci1 =(0.95/2)*(np.std(mat['sens_monosite_h0'][alpha_idx,:,:])/np.sqrt(n))
    ci2 =(0.95/2)*(np.std(mat['sens_multisite_h0_dummy'][alpha_idx,:,:])/np.sqrt(n))
    ci3 =(0.95/2)*(np.std(mat['sens_multisite_h0'][alpha_idx,:,:])/np.sqrt(n))
    print ci1,ci2,ci3
    plt.legend(['Expected'],bbox_to_anchor=(0., 1.05, 1., .002),handletextpad=0.01, loc=3, ncol=3, mode="expand",borderaxespad=0.,frameon=False)

    ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey',
               alpha=0.5)
    ax.xaxis.grid(False)
    plt.title('H0')
    if alpha_idx==0:
        plt.ylim(0,0.01)
    if alpha_idx==1:
        plt.ylim(0,0.1)
    if alpha_idx==2:
        plt.ylim(0,0.5)
        
    plt.tight_layout()
    plt.savefig(path_out,dpi=150)
    #plt.plot(-ci+0.05*np.ones((6,1)),'k:')
    #plt.plot(ci+0.05*np.ones((6,1)),'k:')



    











In [93]:



    


mat['sens_multisite_dummyvar'].shape



    


















    
Out[93]:








(3, 201, 11)

















In [97]:



    


0.05/np.square(201*11)



    


















    
Out[97]:








1.0228042387462385e-08

















In [79]:



    


mat = load_file(main_path + 'figures/default_sim/samps_s120_15pct/results_simu_power.mat')
out_name = path_out + 'h0_s120_15pct_alpha001.png'
plot_boxplot_h0(mat,out_name,120,15)



    


















    






1.02540949223e-05 1.03451627198e-05 9.98642727507e-05










    
























In [80]:



    


mat = load_file(main_path + 'figures/default_sim/samps_s40_50pct/results_simu_power.mat')
out_name = path_out + 'h0_s40_50pct_alpha001.png'
plot_boxplot_h0(mat,out_name,40,50)



    


















    






1.00000294878e-05 9.8982927624e-06 8.75874267523e-06










    
























In [81]:



    


mat = load_file(main_path + 'figures/default_sim/samps_s80_50pct/results_simu_power.mat')
out_name = path_out + 'h0_s80_50pct_alpha001.png'
plot_boxplot_h0(mat,out_name,80,50)



    


















    






9.97260012162e-06 1.01553704836e-05 8.58882618759e-06










    
























In [82]:



    


mat = load_file(main_path + 'figures/default_sim/samps_s120_50pct/results_simu_power.mat')
out_name = path_out + 'h0_s120_50pct_alpha001.png'
plot_boxplot_h0(mat,out_name,120,50)



    


















    






9.73222870619e-06 1.03629748762e-05 8.52935124479e-06










    
























In [83]:



    


mat = load_file(main_path + 'figures/default_sim/samps_s120_30pct/results_simu_power.mat')
out_name = path_out + 'h0_s120_30pct_alpha001.png'
plot_boxplot_h0(mat,out_name,120,30)



    


















    






9.91909230459e-06 1.00561169211e-05 1.4904612681e-05










    
























In [84]:



    


mat = load_file(main_path + 'figures/default_sim/samps_s120_15pct/results_simu_power.mat')
out_name = path_out + 'h0_s120_15pct_alpha001.png'
plot_boxplot_h0(mat,out_name,120,15)



    


















    






1.02540949223e-05 1.03451627198e-05 9.98642727507e-05










    
























In [92]:



    


mat = load_file(main_path + 'figures/default_sim/samps_s120_50pct/results_simu_power.mat')
out_name = path_out + 'h0_s120_50pct_alpha05.png'
plot_boxplot_h0(mat,out_name,120,50,2)



    


















    






6.89082522015e-05 7.03282606469e-05 9.84011621188e-05










    
























In [91]:



    


mat = load_file(main_path + 'figures/default_sim/samps_s120_50pct/results_simu_power.mat')
out_name = path_out + 'h0_s120_50pct_alpha01.png'
plot_boxplot_h0(mat,out_name,120,50,1)



    


















    






3.15937083592e-05 3.08026197087e-05 3.50809769671e-05










    
























In [132]:



    


main_path



    


















    
Out[132]:








'../../../Dropbox/McGill-publication/Papers/Multisite/'

















In [66]:



    


plt.figure(figsize=(4,4))
plt.boxplot([mat['sens_monosite_h0'][2,:,:],mat['sens_multisite_h0_dummy'][2,:,:],mat['sens_multisite_h0'][2,:,:]],notch=1,labels=('Monosite','Multisite dummy','Multisite'))

n = mat['sens_monosite_h0'][2,:,:].reshape(-1,1).shape[0]
ci =(0.95/2)*(np.std(mat['sens_monosite_h0'][2,:,:])/np.sqrt(n))
plt.plot(0.05*np.ones((6,1)),'k--')
#plt.plot(-ci+0.05*np.ones((6,1)),'k:')
#plt.plot(ci+0.05*np.ones((6,1)),'k:')



    


















    
Out[66]:








[<matplotlib.lines.Line2D at 0x7fbf0432f750>]










    
























In [ ]:



    


np.sqrt



    











In [34]:



    


mat['sens_monosite_h0'][2,:,:].reshape(-1,1).shape



    


















    
Out[34]:








(2211, 1)

alpha at 0.01 80% detection power





In [27]:



    


path_tmp = main_path + 'figures/default_sim/samps_s40_50pct/'
m1,s1 = geteffectsize(path_tmp,'results_simu_power.mat',True,1)
path_tmp = main_path + 'figures/default_sim/samps_s80_50pct/'
m2,s2 = geteffectsize(path_tmp,'results_simu_power.mat',True,1) 
path_tmp = main_path + 'figures/default_sim/samps_s120_50pct/'
m3,s3 = geteffectsize(path_tmp,'results_simu_power.mat',True,1) 

path_tmp = main_path + 'figures/default_sim/samps_s40_50pct/'
mm1,sm1 = geteffectsize(path_tmp,'results_simu_power.mat',False,1) 
path_tmp = main_path + 'figures/default_sim/samps_s80_50pct/'
mm2,sm2 = geteffectsize(path_tmp,'results_simu_power.mat',False,1) 
path_tmp = main_path + 'figures/default_sim/samps_s120_50pct/'
mm3,sm3 = geteffectsize(path_tmp,'results_simu_power.mat',False,1)

yref=[]
x = np.linspace(40, 120, 3)
y = np.array([m1,m2,m3])
y2 = np.array([mm1,mm2,mm3])

f1 = interp1d(x, y, kind='linear')
f2 = interp1d(x, y2, kind='linear')

xnew = np.linspace(40, 120, 40)

for i in xnew:
    yref.append((pow.TTestIndPower().solve_power(nobs1=i/2, power=0.8, ratio=1, alpha=0.01)))
yref  = np.array(yref)

plt.figure(figsize=(6,6))
sns.set_context('paper',font_scale=1.8)
plt.plot(xnew,yref,'k--',lw=2)
#plt.scatter(x,y,c="r",s=50)
#plt.scatter(x,y2,c="b",s=50)
plt.errorbar(x, y, yerr=[s1,s2,s3], fmt='rD')
plt.errorbar(x, y2, yerr=[sm1,sm2,sm3], fmt='bD')

plt.legend(['parametric','monosite','multisite (7 sites)'],bbox_to_anchor=(0., 1.05, 1., .002),handletextpad=0.1, loc=3, ncol=3, mode="expand",borderaxespad=0.,frameon=False)
plt.xlabel('Total # subjects')
plt.ylabel('Effect size (Cohen\'s ' + r'$d$)')
plt.grid(True)
axes = plt.gca()
axes.set_xlim([20,140])
axes.set_ylim([0,1.8])
plt.title(r'False-positive rate $\alpha=0.01$')

plt.savefig(path_out + 'samplesize_effectsize_pow80_alpha01.png',dpi=150)

alpha at 0.05 80% detection power





In [28]:



    


path_tmp = main_path + 'figures/default_sim/samps_s40_50pct/'
m1,s1 = geteffectsize(path_tmp,'results_simu_power.mat',True,2)
path_tmp = main_path + 'figures/default_sim/samps_s80_50pct/'
m2,s2 = geteffectsize(path_tmp,'results_simu_power.mat',True,2) 
path_tmp = main_path + 'figures/default_sim/samps_s120_50pct/'
m3,s3 = geteffectsize(path_tmp,'results_simu_power.mat',True,2) 

path_tmp = main_path + 'figures/default_sim/samps_s40_50pct/'
mm1,sm1 = geteffectsize(path_tmp,'results_simu_power.mat',False,2) 
path_tmp = main_path + 'figures/default_sim/samps_s80_50pct/'
mm2,sm2 = geteffectsize(path_tmp,'results_simu_power.mat',False,2) 
path_tmp = main_path + 'figures/default_sim/samps_s120_50pct/'
mm3,sm3 = geteffectsize(path_tmp,'results_simu_power.mat',False,2)

yref=[]
x = np.linspace(40, 120, 3)
y = np.array([m1,m2,m3])
y2 = np.array([mm1,mm2,mm3])

f1 = interp1d(x, y, kind='linear')
f2 = interp1d(x, y2, kind='linear')

xnew = np.linspace(40, 120, 40)

for i in xnew:
    yref.append((pow.TTestIndPower().solve_power(nobs1=i/2, power=0.8, ratio=1, alpha=0.05)))
yref  = np.array(yref)

plt.figure(figsize=(6,6))
sns.set_context('paper',font_scale=1.8)
plt.plot(xnew,yref,'k--',lw=2)
#plt.scatter(x,y,c="r",s=50)
#plt.scatter(x,y2,c="b",s=50)
plt.errorbar(x, y, yerr=[s1,s2,s3], fmt='rD')
plt.errorbar(x, y2, yerr=[sm1,sm2,sm3], fmt='bD')

plt.legend(['parametric','monosite','multisite (7 sites)'],bbox_to_anchor=(0., 1.05, 1., .002),handletextpad=0.1, loc=3, ncol=3, mode="expand",borderaxespad=0.,frameon=False)
plt.xlabel('Total # subjects')
plt.ylabel('Effect size (Cohen\'s ' + r'$d$)')
plt.grid(True)
axes = plt.gca()
axes.set_xlim([20,140])
axes.set_ylim([0,1.8])
plt.title(r'False-positive rate $\alpha=0.05$')
plt.savefig(path_out + 'samplesize_effectsize_pow80_alpha05.png',dpi=150)



    


















    
























In [408]:



    


np.linspace(0, 120, 4)



    


















    
Out[408]:








array([   0.,   40.,   80.,  120.])

Figures for the real data

  • 7 sites of ~20 subjects each, (real data)
  • Balancing of pathology per site is 50%
  • For each figure we vary the sample size from 40, 80 and 120 subjects

conn = 0 path_tmp = main_path + 'figures/samps_s40_50pct/' for conn in range(9, 10): gen_fig(path_tmp,'results_simu_power.mat',path_tmp, 'detection_power.pdf',conn)

path_tmp = main_path + 'figures/samps_s80_50pct/' for conn in range(9, 10): gen_fig(path_tmp,'results_simu_power.mat',path_tmp, 'detection_power.pdf',conn)

path_tmp = main_path + 'figures/samps_s120_50pct/' for conn in range(9, 10): gen_fig(path_tmp,'results_simu_power.mat',path_tmp, 'detection_power.pdf',conn)





In [6]:



    


sns.set_color_codes('muted')
path_tmp = main_path + 'figures/default_sim/samps_s40_50pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s40_50pct.png',efflim=2)


path_tmp = main_path + 'figures/default_sim/samps_s80_50pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s80_50pct.png',efflim=2)


path_tmp = main_path + 'figures/default_sim/samps_s120_50pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s120_50pct.png',efflim=2)


current_palette = sns.color_palette()
sns.palplot(current_palette)

current_palette = sns.color_palette('pastel')
sns.palplot(current_palette)

current_palette = sns.color_palette('muted')
sns.palplot(current_palette)



    


















    






('monosite h0: ', 0.001012211668928087)
('monosite h0: ', 0.00099050203527815484)
('monosite h0: ', 0.00096426956128448668)

Figures for the real data

  • 7 sites of ~20 subjects each, (real data)
  • Balancing of pathology per site is:
    • Figure 1: 15% and 85%
    • Figure 2: 30% and 70%
    • Figure 3: 50% and 50%




In [12]:



    


'''
path_tmp = main_path + 'figures/default_sim/samps_s80_15pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s80_15pct.png')
    
path_tmp = main_path + 'figures/default_sim/samps_s80_30pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s80_30pct.png')
    
path_tmp = main_path + 'figures/default_sim/samps_s80_50pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s80_50pct.png')
'''

path_tmp = main_path + 'figures/default_sim/samps_s120_15pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s120_15pct_14.png',efflim=2)
    
path_tmp = main_path + 'figures/default_sim/samps_s120_30pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s120_30pct_14.png',efflim=2)
    
path_tmp = main_path + 'figures/default_sim/samps_s120_50pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_s120_50pct_14.png',efflim=2)



    


















    






('monosite h0: ', 0.0010194482134780645)
('monosite h0: ', 0.00098914518317503404)
('monosite h0: ', 0.00096426956128448668)

Figures for the real data with site effect in patho

  • 7 sites of ~20 subjects each for a total of 120 subjects, (real data)
  • Effect of the site in relationship with the pathology (2x the effect on site 1 and 1/2 the effect on site 2)
  • Balancing of pathology per site is:
    • Figure 1: 15% and 85%
    • Figure 2: 30% and 70%
    • Figure 3: 50% and 50%




In [19]:



    


path_tmp = main_path + 'figures/sim_site_patho/samps_s80_15pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho_s80_15pct.png',efflim=2)
path_tmp = main_path + 'figures/sim_site_patho/samps_s80_30pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho_s80_30pct.png',efflim=2)
path_tmp = main_path + 'figures/sim_site_patho/samps_s80_50pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho_s80_50pct.png',efflim=2)

path_tmp = main_path + 'figures/sim_site_patho/samps_s120_15pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho_s120_15pct.png',efflim=2)
path_tmp = main_path + 'figures/sim_site_patho/samps_s120_30pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho_s120_30pct.png',efflim=2)
path_tmp = main_path + 'figures/sim_site_patho/samps_s120_50pct/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho_s120_50pct.png',efflim=2)



    


















    






('monosite h0: ', 0.0010104025327905925)
('monosite h0: ', 0.001046132971506106)
('monosite h0: ', 0.00099050203527815484)
('monosite h0: ', 0.0010194482134780645)
('monosite h0: ', 0.00098914518317503404)
('monosite h0: ', 0.00096426956128448668)

Real data simulation 2 sites unbalanced

path_tmp = main_path + 'figures/2sites_realdata/samps_s33_70_30_2site/' gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho2site_s33_70pct.pdf') path_tmp = main_path + 'figures/2sites_realdata/samps_s33_30_70_2site/' gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho2site_s33_30pct.pdf')





In [18]:



    


path_tmp = main_path + 'figures/2sites_realdata/samps_s100_30_70_2site/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho2site_s100_30pct.png',efflim=2,nsites=2)
path_tmp = main_path + 'figures/2sites_realdata/samps_s100_70_30_2site/'
gen_fig_allconn(path_tmp,'results_simu_power.mat',path_out, 'realdata_detect_pow_sitepatho2site_s100_70pct.png',efflim=2,nsites=2)



    


















    






('monosite h0: ', 0.0010004522840343737)
('monosite h0: ', 0.00098190863862505657)

Figures for full simulations

Balanced sites 50/50

  • 2 sites of 50 subjects each,
  • Balancing of pathology per site is 50% and 50%




In [413]:



    


path_tmp = main_path + 'figures/'

gen_fig_full_sim(path_tmp,'simu_bal5050_var0_site0.mat',path_tmp, 'detect_pow_bal5050_var0_site0.png')
gen_fig_full_sim(path_tmp,'simu_bal5050_var0_site05.mat',path_tmp, 'detect_pow_bal5050_var0_site05.png')
gen_fig_full_sim(path_tmp,'simu_bal5050_var2_site0.mat',path_tmp, 'detect_pow_bal5050_var2_site0.png')
gen_fig_full_sim(path_tmp,'simu_bal5050_var2_site05.mat',path_tmp, 'detect_pow_bal5050_var2_site05.png')



    


















    






('monosite h0: ', array([ 0.001]))
('monosite h0: ', array([ 0.00098675]))
('monosite h0: ', array([ 0.00102649]))
('monosite h0: ', array([ 0.00099338]))

Unbalanced site 70/30

  • 2 sites of 50 subjects each,
  • Balancing of pathology per site is 70% and 30%




In [414]:



    


path_tmp = main_path + 'figures/'
gen_fig_full_sim(path_tmp,'simu_bal7030_var0_site0.mat',path_tmp, 'detect_pow_bal7030_var0_site0.png')
gen_fig_full_sim(path_tmp,'simu_bal7030_var0_site05.mat',path_tmp, 'detect_pow_bal7030_var0_site05.png')
gen_fig_full_sim(path_tmp,'simu_bal7030_var2_site0.mat',path_tmp, 'detect_pow_bal7030_var2_site0.png')
gen_fig_full_sim(path_tmp,'simu_bal7030_var2_site05.mat',path_tmp, 'detect_pow_bal7030_var2_site05.png')



    


















    






('monosite h0: ', array([ 0.00102649]))
('monosite h0: ', array([ 0.0011457]))
('monosite h0: ', array([ 0.00102649]))
('monosite h0: ', array([ 0.00095364]))

Unbalanced and sample size difference

  • 2 sites of 20 subjects and 80 subjects,
  • Balancing of pathology per site is 50% and 50%
  • Effect of the site in relationship with the pathology (2x the effect on site 1 and 1/2 the effect on site 2)




In [415]:



    


path_tmp = main_path + 'figures/'
gen_fig_full_sim(path_tmp,'simu_2080bal5050_var0_site0.mat',path_tmp, 'detect_pow_2080bal5050_var0_site0.png')
gen_fig_full_sim(path_tmp,'simu_2080bal5050_var0_site05.mat',path_tmp, 'detect_pow_2080bal5050_var0_site05.png')
gen_fig_full_sim(path_tmp,'simu_2080bal5050_var2_site0.mat',path_tmp, 'detect_pow_2080bal5050_var2_site0.png')
gen_fig_full_sim(path_tmp,'simu_2080bal5050_var2_site05.mat',path_tmp, 'detect_pow_2080bal5050_var2_site05.png')



    


















    






('monosite h0: ', array([ 0.00090066]))
('monosite h0: ', array([ 0.00090066]))
('monosite h0: ', array([ 0.00090728]))
('monosite h0: ', array([ 0.00094702]))
  • 2 sites of 20 subjects and 80 subjects,
  • Balancing of pathology per site is 70% and 30%
  • Effect of the site in relationship with the pathology (2x the effect on site 1 and 1/2 the effect on site 2)




In [416]:



    


path_tmp = main_path + 'figures/'

gen_fig_full_sim(path_tmp,'simu_2080bal7030_var0_site0.mat',path_tmp, 'detect_pow_2080bal7030_var0_site0.png')
gen_fig_full_sim(path_tmp,'simu_2080bal7030_var0_site05.mat',path_tmp, 'detect_pow_2080bal7030_var0_site05.png')
gen_fig_full_sim(path_tmp,'simu_2080bal7030_var2_site0.mat',path_tmp, 'detect_pow_2080bal7030_var2_site0.png')
gen_fig_full_sim(path_tmp,'simu_2080bal7030_var2_site05.mat',path_tmp, 'detect_pow_2080bal7030_var2_site05.png')

gen_fig_full_sim(path_tmp,'simu_8020bal7030_var0_site0.mat',path_tmp, 'detect_pow_8020bal7030_var0_site0.png')
gen_fig_full_sim(path_tmp,'simu_8020bal7030_var0_site05.mat',path_tmp, 'detect_pow_8020bal7030_var0_site05.png')
gen_fig_full_sim(path_tmp,'simu_8020bal7030_var2_site0.mat',path_tmp, 'detect_pow_8020bal7030_var2_site0.png')
gen_fig_full_sim(path_tmp,'simu_8020bal7030_var2_site05.mat',path_tmp, 'detect_pow_8020bal7030_var2_site05.png')



    


















    






('monosite h0: ', array([ 0.00113245]))
('monosite h0: ', array([ 0.00092715]))
('monosite h0: ', array([ 0.00093377]))
('monosite h0: ', array([ 0.00107285]))
('monosite h0: ', array([ 0.0010596]))
('monosite h0: ', array([ 0.00098013]))
('monosite h0: ', array([ 0.00088742]))
('monosite h0: ', array([ 0.00107947]))

Pharma setup

50





In [417]:



    


path_tmp = main_path + 'figures/'

gen_fig_full_sim(path_tmp,'simu_10subj_20sites_rndbal1090_var0_site05.mat',path_tmp, 'detect_pow_10subj_20sites_rndbal1090_var0_site05.png')
gen_fig_full_sim(path_tmp,'simu_10subj_50sites_rndbal1090_var0_site05.mat',path_tmp, 'detect_pow_10subj_50sites_rndbal1090_var0_site05.png')
gen_fig_full_sim(path_tmp,'simu_rndsubj0520_50sites_rndbal1090_var0_site05.mat',path_tmp, 'detect_pow_rndsubj0520_50sites_rndbal1090_var0_site05.png')
gen_fig_full_sim(path_tmp,'simu_rndsubj0215_50sites_rndbal1090_var0_siternd0005.mat',path_tmp, 'detect_pow_rndsubj0215_50sites_rndbal1090_var0_siternd0005.png')
gen_fig_full_sim(path_tmp,'simu_rndsubj0215_50sites_rndbal1090_varrnd02_siternd0005.mat',path_tmp, 'detect_pow_rndsubj0215_50sites_rndbal1090_varrnd02_siternd0005.png')



    


















    






('monosite h0: ', array([ 0.00090728]))
('monosite h0: ', array([ 0.00098013]))
('monosite h0: ', array([ 0.00093377]))
('monosite h0: ', array([ 0.00093377]))
('monosite h0: ', array([ 0.00091391]))










    

















    

















    

















    

















    
























In [418]:



    


aa = load_file(path_tmp + 'simu_10subj_20sites_rndbal1090_var0_site05.mat')



    











In [419]:



    


aa.keys()
aa.get('sens_mono')



    


















    
Out[419]:








array([[ 0.   ],
       [ 0.002],
       [ 0.   ],
       [ 0.003],
       [ 0.   ],
       [ 0.004],
       [ 0.006],
       [ 0.001],
       [ 0.002],
       [ 0.002],
       [ 0.003],
       [ 0.009],
       [ 0.008],
       [ 0.004],
       [ 0.012],
       [ 0.012],
       [ 0.018],
       [ 0.016],
       [ 0.024],
       [ 0.018],
       [ 0.033],
       [ 0.024],
       [ 0.032],
       [ 0.053],
       [ 0.051],
       [ 0.05 ],
       [ 0.059],
       [ 0.074],
       [ 0.092],
       [ 0.093],
       [ 0.112],
       [ 0.129],
       [ 0.132],
       [ 0.165],
       [ 0.195],
       [ 0.177],
       [ 0.218],
       [ 0.234],
       [ 0.253],
       [ 0.264],
       [ 0.295],
       [ 0.338],
       [ 0.35 ],
       [ 0.386],
       [ 0.422],
       [ 0.438],
       [ 0.496],
       [ 0.495],
       [ 0.496],
       [ 0.577],
       [ 0.57 ],
       [ 0.597],
       [ 0.617],
       [ 0.675],
       [ 0.698],
       [ 0.707],
       [ 0.73 ],
       [ 0.733],
       [ 0.747],
       [ 0.784],
       [ 0.802],
       [ 0.864],
       [ 0.866],
       [ 0.868],
       [ 0.89 ],
       [ 0.901],
       [ 0.909],
       [ 0.916],
       [ 0.919],
       [ 0.937],
       [ 0.935],
       [ 0.946],
       [ 0.957],
       [ 0.949],
       [ 0.963],
       [ 0.975],
       [ 0.977],
       [ 0.978],
       [ 0.985],
       [ 0.988],
       [ 0.985],
       [ 0.992],
       [ 0.996],
       [ 0.997],
       [ 0.995],
       [ 0.994],
       [ 0.995],
       [ 0.998],
       [ 0.998],
       [ 0.997],
       [ 0.996],
       [ 0.998],
       [ 0.999],
       [ 0.999],
       [ 1.   ],
       [ 1.   ],
       [ 0.999],
       [ 0.999],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ],
       [ 1.   ]])

















In [420]:



    


#color palette "#9b59b6", "#3498db", "#95a5a6", "#e74c3c", "#34495e", "#2ecc71"
def plotpropbar(patho,ctrl):
    plt.figure(figsize=(2, 2))
    plt.subplot2grid((5,5), (1,1),colspan=4, rowspan=4)
    N = 2
    ind = np.arange(N)/4.    # the x locations for the groups
    width = 0.20       # the width of the bars: can also be len(x) sequence
    p1 = plt.bar(ind, patho,width, color='#3498db')
    p2 = plt.bar(ind, ctrl, width, color='#95a5a6',bottom=patho)

    plt.ylabel('# subjects')
    #plt.title('Scores by group and gender')

    plt.xticks(ind+width/2., ('Site1', 'Site2') )
    plt.yticks(np.arange(0,101,20))
    #plt.legend( (p1[0], p2[0]), ('Patho', 'CTRL') )
    plt.legend((p1[0], p2[0]), ('Patho', 'Control'),bbox_to_anchor=(0., 1.0, 1., .102), loc=3, ncol=1, borderaxespad=0.,frameon=False)

path_tmp = main_path + 'figures/'


# simulation data repartition
patho   = (35, 15)
ctrl = (15, 35)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_5050subj_7030.png')

patho   = (25, 25)
ctrl = (25, 25)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_5050subj_5050.png')

patho   = (40, 10)
ctrl = (40, 10)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_8020subj_5050.png')

patho   = (56, 6)
ctrl    = (24, 14)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_8020subj_7030.png')

patho   = (10, 40)
ctrl = (10, 40)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_2080subj_5050.png')

patho   = (6, 56)
ctrl    = (14, 24)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_2080subj_7030.png')

# real data repartition
patho   = (30, 30)
ctrl    = (30, 30)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_real120subj_5050.png')

patho   = (42, 18)
ctrl    = (18, 42)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_real120subj_7030.png')

patho   = (18, 42)
ctrl    = (42, 18)
plotpropbar(patho,ctrl)
plt.savefig(path_tmp + 'prop_real120subj_3070.png')



    


















    

















    

















    

















    

















    

















    

















    

















    

















    
























In [35]:



    


4990*0.001



    


















    
Out[35]:








4.99

















In [15]:



    


def errorfill_pred(x, y, yerr, color=None, alpha_fill=0.3, ax=None, label=''):
    ax = ax if ax is not None else plt.gca()
    if color is None:
        color = ax._get_lines.color_cycle.next()
    if np.isscalar(yerr) or len(yerr) == len(y):
        ymin = y - yerr
        ymax = y + yerr
    elif len(yerr) == 2:
        ymin, ymax = yerr
    ax.plot(x, y, color=color, label=label, lw=3)
    #ax.fill_between(x, ymax, ymin, color=color, alpha=alpha_fill)
    
def gen_fig_prediction(path_in, file_in, path_out, file_out, conn=0, pvalidx=0, efflim=1.5,nsites=7):
    labels_pct_regions = ['0.1%','1%','5%']
    labels_box = ['H0 multi','H0 Dummy', 'H0 METAL']
    lw = 2
    mat = load_file(path_in + file_in)
    pvals = mat.get('list_p')[0]
    
    fig = plt.figure(figsize=(5,5),dpi=300)
    sns.set_context('paper',font_scale=1.7)
    plt.subplot2grid((5,2), (1,0),colspan=2,rowspan=4)
    ax = fig.gca()
    #ax.set_xticks(np.arange(0,1,0.1))
    #ax.set_yticks(np.arange(0,1.,0.1))
    if efflim==2:
        plt.xlim(0,2)
        effect_size = np.arange(0.0, 2.001, 0.01)
    else:
        plt.xlim(0,1.5)
        effect_size = np.arange(0.0, 1.501, 0.01)
    
    plt.ylabel('Accuracy')
    plt.xlabel('Effect size (Cohen\'s ' + r'$d$)')
    ax.set_ylim((0,1))
    tmpdata = mat.get('sens_monosite')[pvalidx,:,:]
    errorfill_pred(effect_size, tmpdata.mean(axis=1)[:len(effect_size)],tmpdata.std(axis=1)[:len(effect_size)],label="Monosite", color='r')
    tmpdata = mat.get('sens_multisite_nocorr')[pvalidx,:,:]
    errorfill_pred(effect_size, tmpdata.mean(axis=1),tmpdata.std(axis=1),label="Multisite no correction", color='k')
    tmpdata = mat.get('sens_multisite_dummyvar')[pvalidx,:,:]
    errorfill_pred(effect_size, tmpdata.mean(axis=1)[:len(effect_size)],tmpdata.std(axis=1)[:len(effect_size)],label="Multisite ("+str(nsites)+" sites)", color='b')
    #tmpdata = mat.get('sens_multisite_metal')[pvalidx,:,:]
    #errorfill(effect_size, tmpdata.mean(axis=1),tmpdata.std(axis=1),label="Multisite METAL", color='m')
    
    plt.grid(True)
    plt.title(r'Connections affected $\pi_1=$' +  labels_pct_regions[pvalidx])
    
    # fully use the given bounding box.
    plt.legend(bbox_to_anchor=(0., 1.10, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.,frameon=False)
    
    print('monosite h0: ',np.mean(mat.get('sens_monosite_h0')[pvalidx,:,:]))
    print('multisite no corr h0: ',np.mean(mat.get('sens_multisite_h0')[pvalidx,:,:]))
    print('multisite corr h0: ',np.mean(mat.get('sens_multisite_h0_dummy')[pvalidx,:,:]))
    plt.tight_layout()
    plt.savefig(path_out + file_out,dpi=150)



    











In [7]:



    


path_tmp = main_path + 'figures/simulation_prediction/results_simu_power.mat'
mat = load_file(path_tmp)



    


















    






---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-7-fe8ed520f0b9> in <module>()
      1 path_tmp = main_path + 'figures/simulation_prediction/results_simu_power.mat'
----> 2 mat = load_file(path_tmp)

<ipython-input-2-b8847a2d9158> in load_file(path)
     19 
     20 def load_file(path):
---> 21     mat = scipy.io.loadmat(path)
     22     return mat
     23 

//anaconda/lib/python2.7/site-packages/scipy/io/matlab/mio.pyc in loadmat(file_name, mdict, appendmat, **kwargs)
    132     """
    133     variable_names = kwargs.pop('variable_names', None)
--> 134     MR = mat_reader_factory(file_name, appendmat, **kwargs)
    135     matfile_dict = MR.get_variables(variable_names)
    136     if mdict is not None:

//anaconda/lib/python2.7/site-packages/scipy/io/matlab/mio.pyc in mat_reader_factory(file_name, appendmat, **kwargs)
     55        type detected in `filename`.
     56     """
---> 57     byte_stream = _open_file(file_name, appendmat)
     58     mjv, mnv = get_matfile_version(byte_stream)
     59     if mjv == 0:

//anaconda/lib/python2.7/site-packages/scipy/io/matlab/mio.pyc in _open_file(file_like, appendmat)
     21     if isinstance(file_like, string_types):
     22         try:
---> 23             return open(file_like, 'rb')
     24         except IOError as e:
     25             if appendmat and not file_like.endswith('.mat'):

IOError: [Errno 2] No such file or directory: '/Users/christian/Dropbox/McGill-publication/Papers/Multisite/figures/simulation_prediction/results_simu_power.mat'
















In [44]:



    


mat.keys()



    


















    
Out[44]:








['list_effect_size',
 'sens_monosite_h0',
 'sens_multisite_dummyvar',
 'sens_multisite_h0_dummy',
 'sens_multisite_h0',
 'sens_multisite_nocorr',
 'list_p',
 '__header__',
 '__globals__',
 'sens_monosite',
 'sens_monosite_std',
 'sens_multisite_h0_dummy_std',
 'sens_multisite_nocorr_std',
 'sens_multisite_h0_std',
 '__version__',
 'sens_monosite_h0_std',
 'sens_multisite_dummyvar_std']

















In [67]:



    


np.std(mat['sens_monosite'][2,:,:],axis=1)



    


















    
Out[67]:








array([ 0.09071305,  0.09244237,  0.0928979 ,  0.09156976,  0.09288458,
        0.09190983,  0.09523622,  0.0939944 ,  0.09284272,  0.09646136,
        0.09455437,  0.0929668 ,  0.09533122,  0.09630886,  0.09016815,
        0.09058325,  0.0909945 ,  0.0943959 ,  0.09666571,  0.09343484,
        0.09381242,  0.09483632,  0.09205475,  0.09474108,  0.09305456,
        0.09609084,  0.09355191,  0.09181941,  0.0912776 ,  0.0912776 ,
        0.09474573,  0.09200153,  0.09605732,  0.09632978,  0.09442101,
        0.09610793,  0.09552698,  0.09609317,  0.09584106,  0.09865873,
        0.09713875,  0.09720187,  0.09396856,  0.09558994,  0.09741264,
        0.09708926,  0.09628083,  0.09856033,  0.09904132,  0.09617543,
        0.10007925,  0.09802173,  0.0990886 ,  0.09750269,  0.10017061,
        0.10076331,  0.10519797,  0.10331187,  0.10331187,  0.10843153,
        0.10371957,  0.10241662,  0.10664877,  0.1062401 ,  0.10316927,
        0.1070308 ,  0.10630812,  0.10672403,  0.10393885,  0.10409719,
        0.10343833,  0.10736605,  0.10970973,  0.10511973,  0.10572333,
        0.10723495,  0.11370374,  0.10872524,  0.11169994,  0.10789968,
        0.10802437,  0.10936972,  0.11431051,  0.10861757,  0.11086648,
        0.10524945,  0.10756498,  0.11049448,  0.10760158,  0.10501622,
        0.10528462,  0.10367408,  0.10691497,  0.10798077,  0.10557423,
        0.1011282 ,  0.10132802,  0.09953709,  0.09805608,  0.09995047,
        0.09642179,  0.10085745,  0.09774419,  0.09589532,  0.09947053,
        0.09459754,  0.08968248,  0.09300206,  0.08717786,  0.08913066,
        0.08866843,  0.0850555 ,  0.08497812,  0.08230815,  0.08080245,
        0.08038425,  0.08038425,  0.07283697,  0.07412482,  0.07076131,
        0.07071082,  0.06994108,  0.06576464,  0.06349578,  0.06329272,
        0.06208206,  0.063434  ,  0.05720665,  0.05485109,  0.05582131,
        0.05604863,  0.05041962,  0.05353694,  0.05023232,  0.05114354,
        0.04597443,  0.04399467,  0.04227431,  0.04248043,  0.04100174,
        0.03769988,  0.03780507,  0.03546168,  0.0345376 ,  0.03121608,
        0.03043006,  0.02958894,  0.02601488,  0.0275912 ,  0.02386481,
        0.02564054,  0.02417725,  0.02498977,  0.0228899 ,  0.02125537,
        0.02010572,  0.01744229,  0.01726737,  0.01881675,  0.01610172,
        0.01543981,  0.01527265,  0.01345699,  0.01335829,  0.0133409 ,
        0.01424991,  0.01152366,  0.01096283,  0.01127373,  0.01046764,
        0.0109637 ,  0.00838723,  0.01036059,  0.00897019,  0.00803024,
        0.00779447,  0.00648546,  0.00600499,  0.00729962,  0.00661158,
        0.00649764,  0.00642609,  0.00731149,  0.00629828,  0.00725062,
        0.00632092,  0.00581068,  0.00706204,  0.00488529,  0.00462565,
        0.00419609,  0.00474825,  0.00500968,  0.00475738,  0.00513081,
        0.0057049 ,  0.00337287,  0.00373193,  0.00372509,  0.00448125,
        0.00537366])

















In [61]:



    


main_path



    


















    
Out[61]:








'/Users/christian/Dropbox/McGill-publication/Papers/Multisite/'

















In [16]:



    


path_tmp = main_path + 'figures/simulation_prediction/samps_s120_50pct/'
gen_fig_prediction(path_tmp,'results_simu_power.mat',path_tmp, 'prediction_120_01pct.png',pvalidx=0,efflim=2)
path_tmp = main_path + 'figures/simulation_prediction/samps_s80_50pct/'
gen_fig_prediction(path_tmp,'results_simu_power.mat',path_tmp, 'prediction_80_01pct.png',pvalidx=0,efflim=2)



    


















    






('monosite h0: ', 0.49992590618336896)
('multisite no corr h0: ', 0.50003109452736316)
('multisite corr h0: ', 0.50006734186211799)
('monosite h0: ', 0.49996158911325722)
('multisite no corr h0: ', 0.49993444541995907)
('multisite corr h0: ', 0.50001646180860415)










    

















    
























In [17]:



    


path_tmp = main_path + 'figures/simulation_prediction/samps_s120_50pct/'
gen_fig_prediction(path_tmp,'results_simu_power.mat',path_tmp, 'prediction_120_10pct.png',pvalidx=1,efflim=2)
path_tmp = main_path + 'figures/simulation_prediction/samps_s80_50pct/'
gen_fig_prediction(path_tmp,'results_simu_power.mat',path_tmp, 'prediction_80_10pct.png',pvalidx=1,efflim=2)



    


















    






('monosite h0: ', 0.49949040511727083)
('multisite no corr h0: ', 0.50003642501776846)
('multisite corr h0: ', 0.50009186211798173)
('monosite h0: ', 0.49999890254609319)
('multisite no corr h0: ', 0.49998500146327168)
('multisite corr h0: ', 0.50005297044190811)










    

















    
























In [18]:



    


path_tmp = main_path + 'figures/simulation_prediction/samps_s120_50pct/'
gen_fig_prediction(path_tmp,'results_simu_power.mat',path_tmp, 'prediction_120_50pct.png',pvalidx=2,efflim=2)
path_tmp = main_path + 'figures/simulation_prediction/samps_s80_50pct/'
gen_fig_prediction(path_tmp,'results_simu_power.mat',path_tmp, 'prediction_80_50pct.png',pvalidx=2,efflim=2)



    


















    






('monosite h0: ', 0.50009239516702197)
('multisite no corr h0: ', 0.50002949538024166)
('multisite corr h0: ', 0.5001011016346838)
('monosite h0: ', 0.49985367281240844)
('multisite no corr h0: ', 0.50015130231196947)
('multisite corr h0: ', 0.50001455955516538)










    

















    
























In [ ]:

Content source: SIMEXP/Projects

Similar notebooks:

notebook.community | gallery | about