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
|
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
import numpy as np
from nipy.neurospin.glm import glm
from time import time
def display(t, title=''):
import pylab
pylab.figure()
pylab.imshow(t)
pylab.xlabel('Time')
pylab.ylabel('Electrode')
pylab.title(title)
pylab.show()
# Load data
Y = np.load('bb.npy')
nsubjects = Y.shape[2]
nconditions = Y.shape[3]
ndata = nsubjects*nconditions
conditions = np.asarray(range(nconditions))*50 + 50
saturation = np.array([0,0,1,1,1,1])
# Regressors
baseline = np.ones(ndata)
conditions = np.tile(conditions, nsubjects)
saturation = np.tile(saturation, nsubjects)
subject_factor = np.repeat(np.asarray(range(nsubjects)),6)
## Models
#X = np.asarray([baseline, conditions, saturation]).T
#X = np.asarray([conditions, saturation, subject_factor]).T
#formula='y~1+x1+x2+(1|x3)+(x1|x3)+(x2|x3)'
#contrasts = ([0,1,0], [0,0,1])
X = np.asarray([conditions, subject_factor]).T
formula = 'y~x1+(1|x2)'
contrasts = ([1,0], )
# Test: reduce data
y = Y.reshape([Y.shape[0],Y.shape[1],ndata])
y = y[0:3,0:3,:]
# Standard t-stat
print('Starting fitting...')
tic = time()
m = glm(y, X, axis=2, formula=formula, model='mfx')
dt = time()-tic
print(' duration = %d sec' % dt)
m.save('dump')
# Linear contrast
for con in contrasts:
c = m.contrast(con)
t = c.stat()
display(t, title='Linear contrast')
|
