Example neuropower methodology



In [1]:

    
from __future__ import division
from neuropower import cluster
from neuropower import BUM
from neuropower import neuropowermodels
from neuropower import peakdistribution
import nibabel as nib
import nilearn.plotting
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import seaborn as sns
import scipy
import brewer2mpl
import pandas as pd
colors = brewer2mpl.get_map('Set2', 'qualitative', 5).mpl_colors
twocol = brewer2mpl.get_map('Paired', 'qualitative', 5).mpl_colors

%matplotlib inline

load data and extract peaks



In [2]:

    
map_display=nilearn.plotting.plot_stat_map("Zstat1.nii",threshold=2,title='Z-statistic image')



In [3]:

    
SPM = nib.load("Zstat1.nii").get_data()
maskfile = "Mask.nii"
MASK = nib.load(maskfile).get_data()

exc = 2



In [4]:

    
peaks = cluster.PeakTable(SPM,exc,MASK)
pvalues = np.exp(-exc*(np.array(peaks.peak)-exc))
peaks['pval'] = [max(10**(-6),p) for p in pvalues] 
peaks[1:10]

estimate $\pi_1$



In [5]:

    
bum = BUM.EstimatePi1(peaks['pval'],starts=10,seed=1)
bum

sns.set_style("white")
sns.set_palette(sns.color_palette("Paired"))
ax = sns.distplot(peaks['pval'],kde=False,norm_hist=True,bins=10)
ax.set(xlabel="P-values",ylabel="Frequency",title="Histogram of peak p-values")









    Out[5]:





[<matplotlib.text.Text at 0x7f2603bd4588>,
 <matplotlib.text.Text at 0x7f2603b3b550>,
 <matplotlib.text.Text at 0x7f26036ae588>]



In [6]:

    
xn_p = np.arange(0,1,0.01)
alt_p = bum['pi1']*scipy.stats.beta.pdf(xn_p,bum['a'],1)+1-bum['pi1']
nul_p = [1-bum['pi1']]*len(xn_p)



In [7]:

    
ax = sns.distplot(peaks['pval'],kde=False,norm_hist=True,bins=10)
ax.set(xlabel="P-values",ylabel="Frequency",title="Histogram of peak p-values with estimated mixture distribution \n $\pi_0$="+str(1-bum['pi1']))
plt.plot(xn_p,alt_p,lw=3,color=twocol[1])
plt.plot(xn_p,nul_p,lw=3,color=twocol[0])









    Out[7]:





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

Fit model



In [8]:

    
plt.figure(figsize=(4.5, 4))
sns.set_palette(sns.color_palette("Paired"))
ax = sns.distplot(peaks['peak'],kde=False,norm_hist=True,bins=10)
ax.set(xlabel="peak z-values",ylabel="Frequency",title="Histogram of peak z-values (exc = 2.0)")









    Out[8]:





[<matplotlib.text.Text at 0x7f260343ff98>,
 <matplotlib.text.Text at 0x7f2603418a20>,
 <matplotlib.text.Text at 0x7f2603457be0>]



In [9]:

    
modelfit = neuropowermodels.modelfit(peaks.peak,bum['pi1'],exc=exc,starts=10,method="RFT")
modelfit









    Out[9]:





{'maxloglikelihood': 95.73792312423835,
 'mu': 3.8078812608966661,
 'sigma': 1.1722404484182314}



In [10]:

    
xn_t = np.arange(exc,10,0.01)
alt_t = bum['pi1']*neuropowermodels.altPDF(xn_t,modelfit['mu'],modelfit['sigma'],exc=exc)
nul_t = [1-bum['pi1']]*neuropowermodels.nulPDF(xn_t,exc=exc)
mix_t = neuropowermodels.mixPDF(xn_t,pi1=bum['pi1'],mu=modelfit['mu'],sigma=modelfit['sigma'],exc=exc)



In [11]:

    
plt.figure(figsize=(4.5, 4))
ax = sns.distplot(peaks['peak'],kde=False,norm_hist=True,bins=10)
ax.set(xlabel="peak z-values",ylabel="Frequency",title="Histogram of peak z-values with theoretical nul distibution",xlim=[exc,6])
plt.plot(xn_t,nul_t)









    Out[11]:





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



In [12]:

    
plt.figure(figsize=(4.5, 4))
ax = sns.distplot(peaks['peak'],kde=False,norm_hist=True,bins=10)
ax.set(xlabel="peak z-values",ylabel="Frequency",title="Histogram of peak z-values with fitted distribution",xlim=[exc,6])
plt.plot(xn_t,nul_t,color=twocol[0])
plt.plot(xn_t,alt_t,color=twocol[0])
plt.plot(xn_t,mix_t,color=twocol[1])









    Out[12]:





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



In [13]:

    
thresholds = neuropowermodels.threshold(peaks.peak,peaks.pval,FWHM=[8,8,8],voxsize=[3,3,3],nvox=np.sum(MASK),alpha=0.05,exc=exc)
effect_cohen = modelfit['mu']/np.sqrt(13)
power_predicted = []
newsubs = range(10,71)
for s in newsubs:
    projected_effect = effect_cohen*np.sqrt(s)
    powerpred =  {k:1-neuropowermodels.altCDF(v,projected_effect,modelfit['sigma'],exc=exc,method="RFT") for k,v in thresholds.items() if v!='nan'}
    power_predicted.append(powerpred)

power_predicted_df = pd.DataFrame(power_predicted)
thresholds









    Out[13]:





{'BF': 5.6959999999995929,
 'BH': 4.3329999999997426,
 'RFT': 5.5689999999996065,
 'UN': 3.497999999999835}



In [14]:

    
names = ["UN","BF","RFT","BH"]
plt.figure(figsize=(4.5, 4))
ax = sns.distplot(peaks['peak'],kde=False,norm_hist=True,bins=10)
ax.set(xlabel="peak z-values",ylabel="Frequency",title="Histogram of peak z-values with theoretical nul distibution",xlim=[exc,6])
plt.plot(xn_t,nul_t)
k = 0
for method in names:
    k+=1
    plt.axvline(thresholds[method],color=colors[k])



In [15]:

    
plt.figure(figsize=(4.5, 4))
names = ["UN","BF","RFT","BH"]
for col in range(4):
    plt.plot(newsubs,power_predicted_df[names[col]],color=colors[col+1])
plt.legend(loc='lower right')
plt.xlabel("subjects")
plt.ylabel("power")
plt.title("Power curves")









    Out[15]:





<matplotlib.text.Text at 0x7f26031a6400>



In [16]:

    
power_predicted_df['newsubs'] = newsubs
power_predicted_df[0:10]



In [17]:

    
ind = np.min([i for i,elem in enumerate(power_predicted_df['UN']) if elem>0.8])
sub = power_predicted_df['newsubs'][ind]



In [18]:

    
sub









    Out[18]:





18

x	y	z	peak	pval
14.0	53.0	49.0	2.209951	0.657111
16.0	23.0	25.0	4.187122	0.012598
20.0	26.0	29.0	3.869121	0.023796
21.0	26.0	26.0	3.710561	0.032676
21.0	41.0	36.0	2.499141	0.368512
22.0	17.0	18.0	2.848436	0.183256
22.0	40.0	52.0	2.831554	0.189549
25.0	9.0	26.0	5.025146	0.002357
25.0	73.0	53.0	2.173482	0.706831

	BF	BH	RFT	UN	newsubs
0	0.025431	0.227152	0.032748	0.510960	10
1	0.034078	0.265969	0.043307	0.557306	11
2	0.044598	0.306598	0.055970	0.601723	12
3	0.057138	0.348527	0.070861	0.643846	13
4	0.071811	0.391230	0.088059	0.683394	14
5	0.088688	0.434182	0.107591	0.720169	15
6	0.107794	0.476885	0.129432	0.754052	16
7	0.129105	0.518877	0.153507	0.784999	17
8	0.152553	0.559748	0.179691	0.813031	18
9	0.178021	0.599143	0.207812	0.838222	19