<jupyter_start><jupyter_text>Figure: How large should effect sizes be in neuroimaging to have sufficient power?  Specification of alternativeIn 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. Strategy1. Compute the voxelwise threshold for the specified smoothness and volume     * _FweThres = 5.12_2. Compute the required power, so that the omnibus power is 80%     * _PowerPerRegion = 0.928_3. How large should the maximum statistic in a (small) region be to exceed the voxelwise threshold with 0.928 power?     * _muMax = 4.47_5. How does this voxel statistic translate to Cohen's D for a given sample size?    * _See Figure_<jupyter_code>% 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<jupyter_output><empty_output><jupyter_text>1. What is the voxelwise threshold?<jupyter_code># 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))<jupyter_output>ReselSize: 27
Volume: 228483
ReselCount: 8462.33333333
------------
FWE voxelwise GRF threshold: 5.123062<jupyter_text>2. Definition of alternative  Detect 1 regionWe define a 'success' as a situation in which the maximum in the active field exceedsthe threshold.<jupyter_code>Power = 0.8<jupyter_output><empty_output><jupyter_text>Detect 3 regionsWhat if there are three regions and we want to detect all three.What should the power be in each region to have omnibus 80% power (find all 3)<jupyter_code>regions = 3
PowerPerRegion = Power**(1/3)
print("The power in each region should be: "+str(PowerPerRegion))<jupyter_output>The power in each region should be: 0.928317766723<jupyter_text>Detect 1 larger regionThere 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](https://github.com/jokedurnez/neuropower_new_ideas/blob/master/peakdistribution/chengschwartzman_thresholdfree_distribution_simulation.ipynb) for the simulated results.<jupyter_code>PowerLarge = 1-(1-Power)**(1/3)
print("If we expect 3 peaks, then the power per peak should be: "+str(PowerLarge))<jupyter_output>If we expect 3 peaks, then the power per peak should be: 0.415196452357<jupyter_text>3. How large statistic in a field be to exceed the threshold with power 0.80?We quantify this by computing the expected local maximum in the field (which is a null field elevated by value D).We use the distribution of local maxima of Cheng&Schwartzman to compute the power/effect size.<jupyter_code>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))<jupyter_output>The power is sufficient in larger fields if mu equals: 3.17
The power is sufficient for one region if mu equals: 4.0
The power is sufficient for all regions if mu equals: 4.47<jupyter_text>5. From the required voxel statistic to Cohen's D for a given sample size<jupyter_code># 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)
# 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'])
# add jitter for figure
stdev = 0.01*(max(Data.year)-min(Data.year))
Data['year_jitter'] = Data.year+np.random.randn(len(Data))*stdev
# 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'
                       })
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])
    

Medians
# add logscale

Medians.MdSSLog = [np.log(x) for x in Medians.MdSS]
Data.nLog = [np.log(x) for x in Data.n]<jupyter_output><empty_output><jupyter_text>The figure per List (Tal or David)<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>The figure for different Alternatives<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Figure with logscale and only 1 condition<jupyter_code>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")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15035177625293147, 'fraction_numerical': 0.03196386693303222, 'mean_word_length': 5.698080279232112, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 40, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '76426', 'license': 'MIT', 'n_tokens_mistral': 4415, 'n_tokens_neox': 4107, 'n_words': 847, 'path': 'Fig_power/.ipynb_checkpoints/RequiredEffectSizeForSufficientPower-checkpoint.ipynb', 'repo_name': 'poldracklab/power', 'size': '11513'}