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
|
#!/usr/bin/env python3
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
__doc__ = \
"""
This example shows how to get variance and beta estimated from a nipy GLM.
More specifically:
1. A sequence of fMRI volumes are loaded.
2. A design matrix describing all the effects related to the data is computed.
3. A GLM is applied to the dataset, effect and variance images are produced.
Note that this corresponds to a single run.
Needs matplotlib
Author : Bertrand Thirion, 2010--2012
"""
print(__doc__)
from os import getcwd, mkdir, path
import numpy as np
try:
import matplotlib.pyplot as plt
except ImportError:
raise RuntimeError("This script needs the matplotlib library")
# Local import
from get_data_light import DATA_DIR, get_first_level_dataset
from nibabel import Nifti1Image, save
from nipy.labs.viz import cm, plot_map
from nipy.modalities.fmri.design_matrix import make_dmtx
from nipy.modalities.fmri.experimental_paradigm import load_paradigm_from_csv_file
from nipy.modalities.fmri.glm import FMRILinearModel
#######################################
# Data and analysis parameters
#######################################
# volume mask
# This dataset is large
get_first_level_dataset()
data_path = path.join(DATA_DIR, 's12069_swaloc1_corr.nii.gz')
paradigm_file = path.join(DATA_DIR, 'localizer_paradigm.csv')
# timing
n_scans = 128
tr = 2.4
# paradigm
frametimes = np.linspace(0, (n_scans - 1) * tr, n_scans)
# confounds
hrf_model = 'canonical'
drift_model = "cosine"
hfcut = 128
# write directory
write_dir = path.join(getcwd(), 'results')
if not path.exists(write_dir):
mkdir(write_dir)
print(f'Computation will be performed in directory: {write_dir}')
########################################
# Design matrix
########################################
print('Loading design matrix...')
# the example example.labs.write_paradigm_file shows how to create this file
paradigm = load_paradigm_from_csv_file(paradigm_file)['0']
design_matrix = make_dmtx(frametimes, paradigm, hrf_model=hrf_model,
drift_model=drift_model, hfcut=hfcut)
ax = design_matrix.show()
ax.set_position([.05, .25, .9, .65])
ax.set_title('Design matrix')
plt.savefig(path.join(write_dir, 'design_matrix.png'))
dim = design_matrix.matrix.shape[1]
########################################
# Perform a GLM analysis
########################################
print('Fitting a GLM (this takes time)...')
fmri_glm = FMRILinearModel(data_path, design_matrix.matrix,
mask='compute')
fmri_glm.fit(do_scaling=True, model='ar1')
########################################
# Output beta and variance images
########################################
beta_hat = fmri_glm.glms[0].get_beta() # Least-squares estimates of the beta
variance_hat = fmri_glm.glms[0].get_mse() # Estimates of the variance
mask = fmri_glm.mask.get_fdata() > 0
# output beta images
beta_map = np.tile(mask.astype(np.float64)[..., np.newaxis], dim)
beta_map[mask] = beta_hat.T
beta_image = Nifti1Image(beta_map, fmri_glm.affine)
beta_image.get_header()['descrip'] = (
'Parameter estimates of the localizer dataset')
save(beta_image, path.join(write_dir, 'beta.nii'))
print(f"Beta image witten in {write_dir}")
variance_map = mask.astype(np.float64)
variance_map[mask] = variance_hat
# Create a snapshots of the variance image contrasts
vmax = np.log(variance_hat.max())
plot_map(np.log(variance_map + .1),
fmri_glm.affine,
cmap=cm.hot_black_bone,
vmin=np.log(0.1),
vmax=vmax,
anat=None,
threshold=.1, alpha=.9)
plt.show()
|
