In a brain map in an MNI template, with smoothness of 3 times the voxelsize, there are three similar active regions with voxelwise effect size D. The (spatial) size of the regions is relatively small (<200 voxels). We want to know how large D should be in order to have 80% power to detect all three regions using voxelwise FWE thresholding using Random Field Theory. Detect the region means that the maximum in the activated area exceeds the significance threshold.
In [1]:
% matplotlib inline
from __future__ import division
import os
import nibabel as nib
import numpy as np
from neuropower import peakdistribution
import scipy.integrate as integrate
import pandas as pd
import matplotlib.pyplot as plt
import palettable.colorbrewer as cb
In [2]:
# From smoothness + mask to ReselCount
FWHM = 3
ReselSize = FWHM**3
MNI_mask = nib.load("MNI152_T1_2mm_brain_mask.nii.gz").get_data()
Volume = np.sum(MNI_mask)
ReselCount = Volume/ReselSize
print("ReselSize: "+str(ReselSize))
print("Volume: "+str(Volume))
print("ReselCount: "+str(ReselCount))
print("------------")
# From ReselCount to FWE treshold
FweThres_cmd = 'ptoz 0.05 -g %s' %ReselCount
FweThres = os.popen(FweThres_cmd).read()
print("FWE voxelwise GRF threshold: "+str(FweThres))
In [3]:
Power = 0.8
In [4]:
regions = 3
PowerPerRegion = Power**(1/3)
print("The power in each region should be: "+str(PowerPerRegion))
There is about 1 peak in 200 voxels (based on simulations) Before, we made the assumption that the activated field is so small that there is only one local maximum. We expect there to be only one local maximum if the activated region is smaller than 200 voxels (based on simulations), which seems reasonable. What if the activated regions is larger and we expect 3 peaks to appear in the field?
We treat the peaks as independent, which is reasonable, see this notebook for the simulated results.
In [5]:
PowerLarge = 1-(1-Power)**(1/3)
print("If we expect 3 peaks, then the power per peak should be: "+str(PowerLarge))
In [6]:
muRange = np.arange(1.8,5,0.01)
muLarge = []
muSingle = []
for muMax in muRange:
# what is the power to detect a maximum
power = 1-integrate.quad(lambda x:peakdistribution.peakdens3D(x,1),-20,float(FweThres)-muMax)[0]
if power>PowerLarge:
muLarge.append(muMax)
if power>Power:
muSingle.append(muMax)
if power>PowerPerRegion:
muUnion = muMax
break
print("The power is sufficient in larger fields if mu equals: "+str(muLarge[0]))
print("The power is sufficient for one region if mu equals: "+str(muSingle[0]))
print("The power is sufficient for all regions if mu equals: "+str(muUnion))
In [7]:
# Read in data
Data = pd.read_csv("tal_data.txt",sep=" ",header=None,names=['year','n'])
Data['source']='Tal'
David = pd.read_csv("david_data.txt",sep=" ",header=None,names=['year','n'])
David['source']='David'
Data=Data.append(David)
In [8]:
# add detectable effect
Data['deltaUnion']=muUnion/np.sqrt(Data['n'])
Data['deltaSingle']=muSingle[0]/np.sqrt(Data['n'])
Data['deltaLarge']=muLarge[0]/np.sqrt(Data['n'])
In [9]:
# add jitter for figure
stdev = 0.01*(max(Data.year)-min(Data.year))
Data['year_jitter'] = Data.year+np.random.randn(len(Data))*stdev
In [20]:
# Compute medians per year (for smoother)
Medians = pd.DataFrame({'year':
np.arange(start=np.min(Data.year),stop=np.max(Data.year)+1),
'TalMdSS':'nan',
'DavidMdSS':'nan',
'TalMdDSingle':'nan',
'TalMdDUnion':'nan',
'TalMdDLarge':'nan',
'DavidMdDSingle':'nan',
'DavidMdDUnion':'nan',
'DavidMdDLarge':'nan',
'MdSS':'nan',
'DSingle':'nan',
'DUnion':'nan',
'DLarge':'nan'
})
Out[20]:
In [21]:
for yearInd in (range(len(Medians))):
# Compute medians for Tal's data
yearBoolTal = np.array([a and b for a,b in zip(Data.source=="Tal",Data.year==Medians.year[yearInd])])
Medians.TalMdSS[yearInd] = np.median(Data.n[yearBoolTal])
Medians.TalMdDSingle[yearInd] = np.median(Data.deltaSingle[yearBoolTal])
Medians.TalMdDUnion[yearInd] = np.median(Data.deltaUnion[yearBoolTal])
Medians.TalMdDLarge[yearInd] = np.median(Data.deltaLarge[yearBoolTal])
# Compute medians for David's data
yearBoolDavid = np.array([a and b for a,b in zip(Data.source=="David",Data.year==Medians.year[yearInd])])
Medians.DavidMdSS[yearInd] = np.median(Data.n[yearBoolDavid])
Medians.DavidMdDSingle[yearInd] = np.median(Data.deltaSingle[yearBoolDavid])
Medians.DavidMdDUnion[yearInd] = np.median(Data.deltaUnion[yearBoolDavid])
Medians.DavidMdDLarge[yearInd] = np.median(Data.deltaLarge[yearBoolDavid])
# Compute medians for all data
yearBool = np.array(Data.year==Medians.year[yearInd])
Medians.MdSS[yearInd] = np.median(Data.n[yearBool])
Medians.DSingle[yearInd] = np.median(Data.deltaSingle[yearBool])
Medians.DUnion[yearInd] = np.median(Data.deltaUnion[yearBool])
Medians.DLarge[yearInd] = np.median(Data.deltaLarge[yearBool])
In [22]:
Medians
Out[22]:
In [23]:
# add logscale
Medians.MdSSLog = [np.log(x) for x in Medians.MdSS]
Data.nLog = [np.log(x) for x in Data.n]
In [12]:
twocol = cb.qualitative.Paired_12.mpl_colors
fig,axs = plt.subplots(1,2,figsize=(12,5))
fig.subplots_adjust(hspace=.5,wspace=.3)
axs=axs.ravel()
axs[0].plot(Data.year_jitter[Data.source=="Tal"],Data.n[Data.source=="Tal"],"r.",color=twocol[0],alpha=0.5,label="")
axs[0].plot(Data.year_jitter[Data.source=="David"],Data.n[Data.source=="David"],"r.",color=twocol[2],alpha=0.5,label="")
axs[0].plot(Medians.year,Medians.TalMdSS,color=twocol[1],lw=3,label="Neurosynth")
axs[0].plot(Medians.year,Medians.DavidMdSS,color=twocol[3],lw=3,label="David et al.")
axs[0].set_xlim([1993,2016])
axs[0].set_ylim([0,200])
axs[0].set_xlabel("Year")
axs[0].set_ylabel("Median Sample Size")
axs[0].legend(loc="upper left",frameon=False)
axs[1].plot(Data.year_jitter[Data.source=="Tal"],Data.deltaUnion[Data.source=="Tal"],"r.",color=twocol[0],alpha=0.5,label="")
axs[1].plot(Data.year_jitter[Data.source=="David"],Data.deltaUnion[Data.source=="David"],"r.",color=twocol[2],alpha=0.5,label="")
axs[1].plot(Medians.year,Medians.TalMdDUnion,color=twocol[1],lw=3,label="Neurosynth")
axs[1].plot(Medians.year,Medians.DavidMdDUnion,color=twocol[3],lw=3,label="David et al.")
axs[1].set_xlim([1993,2016])
axs[1].set_ylim([0,3])
axs[1].set_xlabel("Year")
axs[1].set_ylabel("Effect Size with 80% power")
axs[1].legend(loc="upper right",frameon=False)
plt.show()
In [17]:
twocol = cb.qualitative.Paired_12.mpl_colors
fig,axs = plt.subplots(1,2,figsize=(12,5))
fig.subplots_adjust(hspace=.5,wspace=.3)
axs=axs.ravel()
axs[0].plot(Data.year_jitter,Data.n,"r.",color=twocol[0],alpha=0.5,label="")
axs[0].plot(Medians.year,Medians.MdSS,color=twocol[1],lw=3,label="Median Sample Size")
axs[0].set_xlim([1993,2016])
axs[0].set_ylim([0,200])
axs[0].set_xlabel("Year")
axs[0].set_ylabel("Sample Size")
axs[0].legend(loc="upper left",frameon=False)
axs[1].plot(Data.year_jitter,Data.deltaUnion,"r.",color=twocol[0],alpha=0.5,label="")
axs[1].plot(Data.year_jitter,Data.deltaSingle,"r.",color=twocol[2],alpha=0.5,label="")
axs[1].plot(Data.year_jitter,Data.deltaLarge,"r.",color=twocol[4],alpha=0.5,label="")
axs[1].plot(Medians.year,Medians.DUnion,color=twocol[1],lw=3,label="Alternative: detect 3 regions")
axs[1].plot(Medians.year,Medians.DSingle,color=twocol[3],lw=3,label="Alternative: detect 1 region")
axs[1].plot(Medians.year,Medians.DLarge,color=twocol[5],lw=3,label="Alternative: detect 1 (larger) region")
axs[1].set_xlim([1993,2016])
axs[1].set_ylim([0,3])
axs[1].set_xlabel("Year")
axs[1].set_ylabel("Effect Size detectable with 80% power")
axs[1].legend(loc="upper right",frameon=False)
#plt.savefig("ReqEffSize.pdf")
Out[17]:
In [27]:
twocol = cb.qualitative.Paired_12.mpl_colors
fig,axs = plt.subplots(1,2,figsize=(12,5))
fig.subplots_adjust(hspace=.5,wspace=.3)
axs=axs.ravel()
axs[0].plot(Data.year_jitter,Data.n,"r.",color=twocol[0],alpha=0.5,label="")
axs[0].plot(Medians.year,Medians.MdSS,color=twocol[1],lw=3,label="log(Median Sample Size)")
axs[0].set_xlim([1993,2016])
axs[0].set_ylim([0,200])
axs[0].set_xlabel("Year")
axs[0].set_ylabel("log(Sample Size)")
axs[0].legend(loc="upper left",frameon=False)
#axs[1].plot(Data.year_jitter,Data.deltaUnion,"r.",color=twocol[0],alpha=0.5,label="")
axs[1].plot(Data.year_jitter,Data.deltaSingle,"r.",color=twocol[2],alpha=0.5,label="")
#axs[1].plot(Data.year_jitter,Data.deltaLarge,"r.",color=twocol[4],alpha=0.5,label="")
#axs[1].plot(Medians.year,Medians.DUnion,color=twocol[1],lw=3,label="Alternative: detect 3 regions")
axs[1].plot(Medians.year,Medians.DSingle,color=twocol[3],lw=3,label="Median detectable effect size")
#axs[1].plot(Medians.year,Medians.DLarge,color=twocol[5],lw=3,label="Alternative: detect 1 (larger) region")
axs[1].set_xlim([1993,2016])
axs[1].set_ylim([0,3])
axs[1].set_xlabel("Year")
axs[1].set_ylabel("Effect Size detectable with 80% power")
axs[1].legend(loc="upper right",frameon=False)
plt.savefig("ReqEffSize.jpg")