1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
Example of script to analyse the reproducibility in group studies
using a bootstrap procedure
This reproduces approximately the work described in
Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses.
Thirion B, Pinel P, Meriaux S, Roche A, Dehaene S, Poline JB.
Neuroimage. 2007 Mar;35(1):105-20.
author: Bertrand Thirion, 2005-2009
"""
import numpy as np
import os.path as op
import cPickle
from nipy.io.imageformats import load, save, Nifti1Image
import tempfile
import get_data_light
from nipy.neurospin.utils.mask import intersect_masks
from nipy.neurospin.utils.reproducibility_measures import \
voxel_reproducibility, cluster_reproducibility, map_reproducibility
print 'This analysis takes a long while, please be patient'
################################################################################
# Set the paths, data, etc.
get_data_light.getIt()
nsubj = 12
subj_id = range(nsubj)
nbeta = 29
data_dir = op.expanduser(op.join('~', '.nipy', 'tests', 'data',
'group_t_images'))
mask_images = [op.join(data_dir,'mask_subj%02d.nii'%n)
for n in range(nsubj)]
stat_images =[ op.join(data_dir,'spmT_%04d_subj_%02d.nii'%(nbeta,n))
for n in range(nsubj)]
contrast_images =[ op.join(data_dir,'con_%04d_subj_%02d.nii'%(nbeta,n))
for n in range(nsubj)]
swd = tempfile.mkdtemp('image')
################################################################################
# Make a group mask
affine = load(mask_images[0]).get_affine()
mask = intersect_masks(mask_images)>0
grp_mask = Nifti1Image(mask, affine)
xyz = np.where(mask)
xyz = np.array(xyz).T
nvox = xyz.shape[0]
################################################################################
# Load the functional images
# Load the betas
Functional = []
VarFunctional = []
tiny = 1.e-15
for s in range(nsubj):
beta = []
varbeta = []
rbeta = load(contrast_images[s])
temp = (rbeta.get_data())[mask]
beta.append(temp)
rbeta = load(stat_images[s])
tempstat = (rbeta.get_data())[mask]
# ugly trick to get the variance
varbeta = temp**2/(tempstat**2+tiny)
varbeta = np.array(varbeta)
beta = np.array(beta)
Functional.append(beta.T)
VarFunctional.append(varbeta.T)
Functional = np.array(Functional)
Functional = np.squeeze(Functional).T
VarFunctional = np.array(VarFunctional)
VarFunctional = np.squeeze(VarFunctional).T
Functional[np.isnan(Functional)] = 0
VarFunctional[np.isnan(VarFunctional)] = 0
################################################################################
# script
ngroups = 10
thresholds = [1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0]
sigma = 6.0
csize = 10
niter = 10
method = 'crfx'
verbose = 0
# apply random sign swaps to the data to test the reproducibility
# under the H0 hypothesis
swap = False
kap = []
clt = []
for threshold in thresholds:
kappa = []
cls = []
kwargs={'threshold':threshold,'csize':csize}
for i in range(niter):
k = voxel_reproducibility(Functional, VarFunctional, grp_mask, ngroups,
method, swap, verbose, **kwargs)
kappa.append(k)
cld = cluster_reproducibility(Functional, VarFunctional, grp_mask, ngroups,
sigma, method, swap,
verbose, **kwargs)
cls.append(cld)
kap.append(np.array(kappa))
clt.append(np.array(cls))
import matplotlib.pylab as mp
mp.figure()
mp.subplot(1,2,1)
mp.boxplot(kap)
mp.title('voxel-level reproducibility')
mp.xticks(range(1,1+len(thresholds)),thresholds)
mp.xlabel('threshold')
mp.subplot(1,2,2)
mp.boxplot(clt)
mp.title('cluster-level reproducibility')
mp.xticks(range(1,1+len(thresholds)),thresholds)
mp.xlabel('threshold')
################################################################################
# create an image
th = 4.0
swap = False
kwargs = {'threshold':th,'csize':csize}
rmap = map_reproducibility(Functional, VarFunctional, grp_mask, ngroups,
method, swap, verbose, **kwargs)
wmap = mask.astype(np.int)
wmap[mask] = rmap
wim = Nifti1Image(wmap, affine)
wim.get_header()['descrip']= 'reproducibility map at threshold %f, \
cluster size %d'%(th,csize)
wname = op.join(swd,'repro.nii')
save(wim, wname)
print('Wrote a reproducibility image in %s'%wname)
|
