A simple second level t-test using nipy¶
Simple random effect model¶
Random effect model in neuroimaging refer to the fact that the variance used for testing the effect is the between subject variance. In other words, we take the effect computed for each subject (the contrast maps) and we enter these in a second level effect model.
# try Bertrand's and Alexis example on Matthew's data :
import numpy as np
import nibabel as nib
from nibabel import Nifti1Image as Image
import nipy
#from nipy.utils import example_data
print nipy.__file__ # see which nipy we are using
Get the data located in ./ds107/analysis/sub0*.img
import os
import os.path as osp
from os.path import join as pjoin, isfile
import zipfile as zpf
#----- SET THE DIRECTORY WHERE YOUR DATA SIT ------#
DATADIR = pjoin(osp.expanduser('~'), 'data', 'ds107')
#----- if DATA ARE ZIP : UNZIP --------------------#
datazip = pjoin(pjoin(DATADIR,'analysis.zip'))
if osp.isfile(datazip):
with zpf.ZipFile(datazip, "r") as z:
z.extractall(DATADIR)
list_subj = range(1,11)
dsubj = [pjoin(DATADIR,'analysis','sub' + "%03d" % i) for i in list_subj]
f_cnames = [pjoin(dsubj[i-1],'copes','_model_id_1cope01.nii.gz') for i in list_subj]
f_vnames = [pjoin(dsubj[i-1],'varcopes','_model_id_1varcope01.nii.gz') for i in list_subj]
#fnames = [os.path.join(scans_dir, 'snn03055dy%i.img' % i) for i in range(1,13)]
%pwd
print [ (f, isfile(f)) for f in f_cnames]
# keep only images that are on disk !
f_cnames = [f for f in f_cnames if isfile(f)]
f_vnames = [f for f in f_vnames if isfile(f)]
c_images = [nib.load(f) for f in f_cnames if isfile(f)]
v_images = [nib.load(f) for f in f_vnames if isfile(f)]
print len(c_images), len(v_images)
Computing the mask of our images - within and between subjects¶
from nipy.labs.mask import compute_mask, compute_mask_files
# first try with the contrast image:
image = c_images[0]
def plot_mask_image(image):
data_img0 = image.get_data()
aff0 = image.get_affine()
hdr0 = image.get_header()
mask_arr = compute_mask(data_img0, None, m=0.2, M=0.9, \
cc=True, opening=2, exclude_zeros=False)
mask_img = nib.Nifti1Image(mask_arr, aff0, hdr0)
f, (a1, a2) = subplots(1,2,figsize=(10,3))
a1.imshow(data_img0[:,:,12])
a1.set_title('my data example')
a2.imshow(mask_arr[:,:,12])
a2.set_title('my mask example')
plot_mask_image(image)
image = v_images[0]
plot_mask_image(image)
the_data = image.get_data()
the_data.shape
#plt.hist(the_data.ravel())
plt.hist(np.log(the_data.ravel() + 1),50)
Usually, you will like to use all subjects (for instance the mean across subjects) to compute the mask
print f_vnames
# you can also compute the mask from the files names (or list of files) directly
mask_arr, mean_img = compute_mask_files(f_vnames, output_filename=None, \
return_mean=True, m=0.2, M=0.9, cc=1, \
exclude_zeros=False, opening=2)
mask_img = nib.Nifti1Image(mask_arr, image.get_affine(), image.get_header())
# print type(mask_arr)
f, (a1, a2) = subplots(1,2,figsize=(10,3))
a1.imshow(mean_img[:,:,12])
a1.set_title('my_example')
a2.imshow(mask_arr[:,:,12])
Construct an "interesting" model:¶
image = c_images[0]
images = c_images
# transformation from voxel no to mm
M = image.get_affine()
# mm to voxel no
iM = np.linalg.inv(M)
# coordinates of voxel of interest in mm (MNI space)
posmm = [-20.0, -42, 34]
# coordinates in voxel space (in homogenous coordinates)
posvox = np.dot(iM, posmm + [1])
# We grab the spatial part of the output. Since we want to use it as an
# index, we need to make it a tuple
posvox = tuple(posvox[:3].astype(int))
print posvox
nimgs = len(images)
vdata = np.zeros(nimgs)
gdata = np.zeros(nimgs)
for i, V in enumerate(images):
d = V.get_data()
vdata[i] = d[posvox]
s = np.std(vdata)/20
x1 = vdata + np.random.randn(vdata.shape[0])*s
print x1
X = np.hstack((vdata[:,newaxis], ones((len(images),1))))
print X
Checking that all files have the same affine :¶
aff0 = v_images[0].get_affine()
shape0 = v_images[0].shape
for i,img in enumerate(v_images):
if i < 2: # print only ..
print img.get_data().shape
print img.get_affine()
print "same as the first one ?: ", np.allclose(aff0, img.get_affine())
#go from a list of 3d to one 4d image:
data_images_4d = nib.concat_images(c_images, check_affines=False)
from nipy.modalities.fmri.glm import FMRILinearModel
#data_images = [nib.load(f) for f in f_cnames]
# FMRILinearModel will
# * load the data in fmri_model.fmri_data
# * load X in fmri_model.design_matrices
# * load the mask or compute it from the runs
fmri_model = FMRILinearModel(data_images_4d, X, mask_img)
# GLM fitting using ordinary least_squares
fmri_model.fit(do_scaling=False, model='ols')
# Here, we were doing the first level analysis, we would generally leave this unspecified,
# and the fit would do a first pass with model="ar1", whiten the data, then a second pass
# with model = 'ols'
# specify and estimate the contrast
contrast_val = np.array(([[1, 0]])) #
print contrast_val
# By default would do a t-test
z_map, = fmri_model.contrast(contrast_val, contrast_type='t', con_id='', output_z=True)
from nipy.utils import templates, DataError
templates
# look at the result
from nipy.labs.viz import plot_map, cm
import matplotlib.pyplot as plt
#anat =
vmax = max(- z_map.get_data().min(), z_map.get_data().max())
vmin = - vmax
plot_map(z_map.get_data(), z_map.get_affine(), #anat=anat,
cut_coords=(-20.0, -42, 34),
cmap=cm.cold_hot,
vmin=vmin,
vmax=vmax,
threshold=1.5,
black_bg=True)
#plt.show()
By looking at the object, we see a number of useful options:
# do something like
fmri_model.contrast??
F contrast can be handled the same way. If you have several lines on your contrast, by defaut this will be treated as an F-test.
# For instance :
contrast_val = np.array(([ [[1, 0],[0, 1]] ])) #
Here are some: ....
- output_z : bool, optional
Return or not the corresponding z-stat image - output_stat : bool, optional
Return or not the base (t/F) stat image - output_effects : bool, optional
Return or not the corresponding effect image - output_variance : bool, optional
Return or not the corresponding variance image
Mixed effect.¶
Mixed effect models refer to the situation when in a model some of the effects are random, and some are fixed. In neuroimaging, the issue is that there is two variances to think of.
- The intra run variance
- The inter subject variance.
The problem:¶
When you want to estimate the variance between subject, you are dealing with noisy observation: the total variance has to be decomposed. We often have a good idea of the first level variance, because in general there are many scans per run, but still there's no direct formula to get the intersubject variance.
For one voxel, let's $Y$ be the observed effects (a vector of $n$-subjects components)
As our $Y$ are measured with a know variance, we have
$$ Y = Z + \epsilon_1 $$where the $Z_i$ represents the "true" measurement for suject $i$, $\epsilon_1 \sim N(0,\Sigma)$, with $\Sigma = diag(\sigma_i^2)$
and at the group level, we have a model between subjects for the $Z$ :
$$ Z = X\beta + \epsilon_2 $$where $\beta$ is the group effect, and $\epsilon_2 \sim N(0,vI)$, $v$ being the group variance (subjects are assumed to be independant).
Overall, we have
$\mathrm{var}(Y) = \Sigma + vI_n $
Most of good packages now include a way of dealing with this. While it is possible to estimate simultaneously $v$, $\Sigma$, and $\beta$ , this requires lenghty computations, and we generally have a good confidence on the $\sigma_i$, such that they are assumed known after a first level estimation, so that only $v$ and $\beta$ remains to be estimated.
here - an explanation of the EM ?¶
Nipy has a couple of routines for this. Let's examine one with a simple example. First create an array of data:
n_sub = 15 # the number of subjects
n_vox = 77 # the number of voxels
# inter subject variance
V2 = np.arange(1,2,1.0/n_sub)
V2 = V2[:,np.newaxis]
vox_variance = np.arange(1.,10.,9.0/n_vox)
vox_variance = vox_variance[np.newaxis,:]
V2 = np.dot(V2, vox_variance)
# within subject variance > within variance
V1 = .5*np.ones(V2.shape)
print np.max(V1)
# check that the variances are positive
if (V1 < 0).any():
raise ValueError('Variance should be positive')
Y = np.random.standard_normal(V1.shape)
Y *= np.sqrt(V2 + V1)
X = np.ones(n_sub)
if X.ndim == 1:
X = X[:, np.newaxis]
# put enough signal such that we should see something significant
beta = 1.
if np.isscalar(beta):
beta = beta * np.ones((X.shape[1], V1.shape[1]))
if beta.ndim == 1:
beta = beta[np.newaxis]
Y += np.dot(X, beta)
print np.max(Y)
from nipy.algorithms.statistics.mixed_effects_stat import mfx_stat
# column to be tested in X
column = 0
t, varh, betah = mfx_stat(Y, V1, X, column, n_iter=5, return_t=True, return_f=False, return_effect=True, \
return_var=True, verbose=False)
print t.shape, betah.shape, varh.shape
print sum(t>2.)
fig = plt.figure()
f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharex=True, figsize=(12,4))
ax1.plot(t)
ax2.plot(betah)
ax3.plot(varh)
Calculate significant voxels and clusters : see Multiple comparison notebook¶