Skip to content

Commit af6cbca

Browse files
matthew-brettmatthew-brett
committed
RF+TST: refactor formula/contrasts to helper func
Move creation of formula and contrasts for event and block designs into a helper function. Add tests for block_design, add creation of dummy block type when no factor given, and fix bug for single block factor.
1 parent fb9bdc1 commit af6cbca

File tree

2 files changed

+177
-127
lines changed

2 files changed

+177
-127
lines changed

‎nipy/modalities/fmri/design.py

Lines changed: 48 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,43 @@ def natural_spline(tvals, knots=None, order=3, intercept=True):
8383
return f.design(tvals, return_float=True)
8484

8585

86+
def _build_formula_contrasts(spec, fields, order):
87+
""" Build formula and contrast in event / block space
88+
89+
Parameters
90+
----------
91+
spec : stuctured array
92+
Structured array containing at least fields listed in `fields`.
93+
fields : sequence of str
94+
Sequence of field names containing names of factors.
95+
order : int
96+
Maximum order of interactions between main effects.
97+
98+
Returns
99+
-------
100+
e_factors : :class:`Formula` instance
101+
Formula for factors given by `fields`
102+
e_contrasts : dict
103+
Dictionary containing contrasts of main effects and interactions
104+
between factors.
105+
"""
106+
e_factors = [Factor(n, np.unique(spec[n])) for n in fields]
107+
e_formula = np.product(e_factors)
108+
e_contrasts = {}
109+
# Add contrasts for factors and factor interactions
110+
max_order = min(len(e_factors), order)
111+
for i in range(1, max_order + 1):
112+
for comb in combinations(zip(fields, e_factors), i):
113+
names = [c[0] for c in comb]
114+
# Collect factors where there is more than one level
115+
fs = [fc.main_effect for fn, fc in comb if len(fc.levels) > 1]
116+
if len(fs) > 0:
117+
e_contrast = np.product(fs).design(spec)
118+
e_contrasts[":".join(names)] = e_contrast
119+
e_contrasts['constant'] = formulae.I.design(spec)
120+
return e_formula, e_contrasts
121+
122+
86123
def event_design(event_spec, t, order=2, hrfs=(glover,),
87124
level_contrasts=False):
88125
""" Create design matrix from event specification `event_spec`
@@ -131,27 +168,13 @@ def event_design(event_spec, t, order=2, hrfs=(glover,),
131168
itertools.cycle([1])),
132169
('time', '_event_'))
133170
fields = ['_event_']
134-
e_factors = [Factor(n, np.unique(event_spec[n])) for n in fields]
135-
e_formula = np.product(e_factors)
136-
137-
# Design and contrasts in event space
171+
e_formula, e_contrasts = _build_formula_contrasts(
172+
event_spec, fields, order)
173+
# Design and contrasts in block space
138174
# TODO: make it so I don't have to call design twice here
139175
# to get both the contrasts and the e_X matrix as a recarray
140176
e_X = e_formula.design(event_spec)
141177
e_dtype = e_formula.dtype
142-
e_contrasts = {}
143-
144-
# Add contrasts for factors and factor interactions
145-
max_order = min(len(e_factors), order)
146-
for i in range(1, max_order + 1):
147-
for comb in combinations(zip(fields, e_factors), i):
148-
names = [c[0] for c in comb]
149-
# Collect factors where there is more than one level
150-
fs = [fc.main_effect for fn, fc in comb if len(fc.levels) > 1]
151-
if len(fs) > 0:
152-
e_contrast = np.product(fs).design(event_spec)
153-
e_contrasts[":".join(names)] = e_contrast
154-
e_contrasts['constant'] = formulae.I.design(event_spec)
155178

156179
# Now construct the design in time space
157180
t_terms = []
@@ -229,40 +252,27 @@ def block_design(block_spec, t, order=2, hrfs=(glover,),
229252
raise ValueError('expecting fields called "start" and "end"')
230253
fields.pop(fields.index('start'))
231254
fields.pop(fields.index('end'))
232-
e_factors = [Factor(n, np.unique(block_spec[n])) for n in fields]
233-
e_formula = np.product(e_factors)
234-
e_contrasts = {}
235-
236-
if len(e_factors) > 1:
237-
for i in range(1, order+1):
238-
for comb in combinations(zip(fields, e_factors), i):
239-
names = [c[0] for c in comb]
240-
fs = [c[1].main_effect for c in comb]
241-
e_contrasts[":".join(names)] = np.product(fs).design(block_spec)
242-
else:
243-
# only one factor, produce the main effect
244-
field = fields[0]
245-
factor = e_factors[0]
246-
e_contrasts[field + '_effect'] = factor.main_effect.design(event_spec)
247-
248-
e_contrasts['constant'] = formulae.I.design(block_spec)
249-
255+
if len(fields) == 0: # No factors specified, make generic block
256+
block_spec = make_recarray(zip(block_spec['start'],
257+
block_spec['end'],
258+
itertools.cycle([1])),
259+
('start', 'end', '_block_'))
260+
fields = ['_block_']
261+
e_formula, e_contrasts = _build_formula_contrasts(
262+
block_spec, fields, order)
250263
# Design and contrasts in block space
251264
# TODO: make it so I don't have to call design twice here
252265
# to get both the contrasts and the e_X matrix as a recarray
253-
254266
e_X = e_formula.design(block_spec)
255267
e_dtype = e_formula.dtype
256268

257269
# Now construct the design in time space
258-
259270
block_times = np.array([(s,e) for s, e in zip(block_spec['start'],
260271
block_spec['end'])])
261272
convolution_interval = (block_times.min() - convolution_padding,
262273
block_times.max() + convolution_padding)
263274

264275
t_terms = []
265-
t_names = []
266276
t_contrasts = {}
267277
for l, h in enumerate(hrfs):
268278
for n in e_dtype.names:

‎nipy/modalities/fmri/tests/test_design.py

Lines changed: 129 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,8 @@
44
import numpy as np
55

66
from ..design import event_design, block_design
7-
from ..utils import events, lambdify_t
7+
from ..utils import (events, lambdify_t, T, convolve_functions,
8+
blocks)
89
from ..hrf import glover
910
from nipy.algorithms.statistics.formula import make_recarray
1011

@@ -26,95 +27,134 @@ def test_event_design():
2627
# Test event design helper function
2728
# An event design with one event type
2829
onsets = np.array([0, 20, 40, 60])
30+
durations = np.array([2, 3, 4, 5])
31+
offsets = onsets + durations
2932
c_fac = np.array([1, 1, 1, 1]) # constant factor (one level)
3033
fac_1 = np.array([0, 1, 0, 1]) # factor 1, two levels
3134
fac_2 = np.array([0, 0, 1, 1]) # factor 2, two levels
3235
t = np.arange(0, 100, 2)
33-
# Event spec with no event factor -> single column design, no contrasts
34-
event_spec_0 = make_recarray(onsets, ('time',))
35-
X_0, contrasts_0 = event_design(event_spec_0, t)
36-
exp_x_0 = lambdify_t(events(onsets, f=glover))(t)
37-
assert_almost_equal(X_0, exp_x_0)
38-
assert_dict_almost_equal(contrasts_0, dict(constant_0=1))
39-
X_0, contrasts_0 = event_design(event_spec_0, t, level_contrasts=True)
40-
assert_almost_equal(X_0, exp_x_0)
41-
assert_dict_almost_equal(contrasts_0, dict(constant_0=1, _event__1_0=1))
42-
# Event spec with single factor, but only one level
43-
event_spec_1c = make_recarray(zip(onsets, c_fac), ('time', 'smt'))
44-
X_1c, contrasts_1c = event_design(event_spec_1c, t)
45-
assert_almost_equal(X_1c, exp_x_0)
46-
assert_dict_almost_equal(contrasts_1c, dict(constant_0=1))
47-
X_1c, contrasts_1c = event_design(event_spec_1c, t, level_contrasts=True)
48-
assert_dict_almost_equal(contrasts_1c, dict(constant_0=1, smt_1_0=1))
49-
# Event spec with single factor, two levels
50-
event_spec_1d = make_recarray(zip(onsets, fac_1), ('time', 'smt'))
51-
exp_x_0 = lambdify_t(events(onsets[fac_1 == 0], f=glover))(t)
52-
exp_x_1 = lambdify_t(events(onsets[fac_1 == 1], f=glover))(t)
53-
X_1d, contrasts_1d = event_design(event_spec_1d, t)
54-
assert_almost_equal(X_1d, np.c_[exp_x_0, exp_x_1])
55-
assert_dict_almost_equal(contrasts_1d,
56-
dict(constant_0=[1, 1], smt_0=[1, -1]))
57-
X_1d, contrasts_1d = event_design(event_spec_1d, t, level_contrasts=True)
58-
assert_dict_almost_equal(contrasts_1d,
59-
dict(constant_0=1,
60-
smt_0=[1, -1], # main effect
61-
smt_0_0=[1, 0], # level 0, hrf 0
62-
smt_1_0=[0, 1])) # level 1, hrf 0
63-
# Event spec with two factors, one with two levels, another with one
64-
event_spec_2dc = make_recarray(zip(onsets, fac_1, c_fac),
65-
('time', 'smt', 'smte'))
66-
X_2dc, contrasts_2dc = event_design(event_spec_2dc, t)
67-
assert_almost_equal(X_2dc, np.c_[exp_x_0, exp_x_1])
68-
assert_dict_almost_equal(contrasts_2dc,
69-
{'constant_0': [1, 1],
70-
'smt_0': [1, -1], # main effect
71-
'smt:smte_0': [1, -1], # interaction
72-
})
73-
X_2dc, contrasts_2dc = event_design(event_spec_2dc, t, level_contrasts=True)
74-
assert_dict_almost_equal(contrasts_2dc,
75-
{'constant_0': [1, 1],
76-
'smt_0': [1, -1], # main effect
77-
'smt:smte_0': [1, -1], # interaction
78-
'smt_0*smte_1_0': [1, 0], # smt 0, smte 0, hrf 0
79-
'smt_1*smte_1_0': [0, 1], # smt 1, smte 0, hrf 0
80-
})
81-
# Event spec with two factors, both with two levels
82-
event_spec_2dd = make_recarray(zip(onsets, fac_1, fac_2),
83-
('time', 'smt', 'smte'))
84-
exp_x_0 = lambdify_t(events(onsets[(fac_1 == 0) & (fac_2 == 0)], f=glover))(t)
85-
exp_x_1 = lambdify_t(events(onsets[(fac_1 == 0) & (fac_2 == 1)], f=glover))(t)
86-
exp_x_2 = lambdify_t(events(onsets[(fac_1 == 1) & (fac_2 == 0)], f=glover))(t)
87-
exp_x_3 = lambdify_t(events(onsets[(fac_1 == 1) & (fac_2 == 1)], f=glover))(t)
88-
X_2dd, contrasts_2dd = event_design(event_spec_2dd, t)
89-
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
90-
exp_cons = {'constant_0': [1, 1, 1, 1],
91-
'smt_0': [1, 1, -1, -1], # main effect fac_1
92-
'smte_0': [1, -1, 1, -1], # main effect fac_2
93-
'smt:smte_0': [1, -1, -1, 1], # interaction
94-
}
95-
assert_dict_almost_equal(contrasts_2dd, exp_cons)
96-
X_2dd, contrasts_2dd = event_design(event_spec_2dd, t, level_contrasts=True)
97-
level_cons = exp_cons.copy()
98-
level_cons.update({
99-
'smt_0*smte_0_0': [1, 0, 0, 0], # smt 0, smte 0, hrf 0
100-
'smt_0*smte_1_0': [0, 1, 0, 0], # smt 0, smte 1, hrf 0
101-
'smt_1*smte_0_0': [0, 0, 1, 0], # smt 1, smte 0, hrf 0
102-
'smt_1*smte_1_0': [0, 0, 0, 1], # smt 1, smte 1, hrf 0
103-
})
104-
assert_dict_almost_equal(contrasts_2dd, level_cons)
105-
# Test max order >> 2, no error
106-
X_2dd, contrasts_2dd = event_design(event_spec_2dd, t, order=100)
107-
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
108-
assert_dict_almost_equal(contrasts_2dd, exp_cons)
109-
# Test max order = 1
110-
X_2dd, contrasts_2dd = event_design(event_spec_2dd, t, order=1)
111-
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
112-
# No interaction
113-
assert_dict_almost_equal(contrasts_2dd,
114-
{'constant_0': [1, 1, 1, 1],
115-
'smt_0': [1, 1, -1, -1], # main effect fac_1
116-
'smte_0': [1, -1, 1, -1], # main effect fac_2
117-
})
118-
# Test field called "time" is necessary
119-
event_spec_1d = make_recarray(zip(onsets, fac_1), ('brighteyes', 'smt'))
120-
assert_raises(ValueError, event_design, event_spec_1d, t)
36+
37+
def mk_ev_spec(factors, names):
38+
names = ('time',) + tuple(names)
39+
if len(factors) == 0:
40+
return make_recarray(onsets, names)
41+
return make_recarray(zip(onsets, *factors), names)
42+
43+
def mk_blk_spec(factors, names):
44+
names = ('start', 'end') + tuple(names)
45+
return make_recarray(zip(onsets, offsets, *factors), names)
46+
47+
def mk_ev_tc(ev_inds):
48+
# Make real time course for given event onsets
49+
return lambdify_t(events(onsets[ev_inds], f=glover))(t)
50+
51+
def mk_blk_tc(ev_inds):
52+
# Make real time course for block onset / offsets
53+
B = blocks(zip(onsets[ev_inds], offsets[ev_inds]))
54+
term = convolve_functions(B, glover(T),
55+
(-5, 70), # step func support
56+
(0, 30.), # conv kernel support
57+
0.02) # dt
58+
return lambdify_t(term)(t)
59+
60+
for d_maker, spec_maker, tc_maker, null_name in (
61+
(event_design, mk_ev_spec, mk_ev_tc, '_event_'),
62+
(block_design, mk_blk_spec, mk_blk_tc, '_block_')):
63+
# Event spec with no event factor -> single column design, no contrasts
64+
spec_0 = spec_maker((), ())
65+
X_0, contrasts_0 = d_maker(spec_0, t)
66+
exp_x_0 = tc_maker(onsets==onsets)
67+
assert_almost_equal(X_0, exp_x_0)
68+
assert_dict_almost_equal(contrasts_0, dict(constant_0=1))
69+
X_0, contrasts_0 = d_maker(spec_0, t, level_contrasts=True)
70+
assert_almost_equal(X_0, exp_x_0)
71+
assert_dict_almost_equal(contrasts_0,
72+
{'constant_0': 1,
73+
null_name + '_1_0': 1})
74+
# Event spec with single factor, but only one level
75+
spec_1c = spec_maker((c_fac,), ('smt',))
76+
X_1c, contrasts_1c = d_maker(spec_1c, t)
77+
assert_almost_equal(X_1c, exp_x_0)
78+
assert_dict_almost_equal(contrasts_1c, dict(constant_0=1))
79+
X_1c, contrasts_1c = d_maker(spec_1c, t, level_contrasts=True)
80+
assert_dict_almost_equal(contrasts_1c, dict(constant_0=1, smt_1_0=1))
81+
# Event spec with single factor, two levels
82+
spec_1d = spec_maker((fac_1,), ('smt',))
83+
exp_x_0 = tc_maker(fac_1 == 0)
84+
exp_x_1 = tc_maker(fac_1 == 1)
85+
X_1d, contrasts_1d = d_maker(spec_1d, t)
86+
assert_almost_equal(X_1d, np.c_[exp_x_0, exp_x_1])
87+
assert_dict_almost_equal(contrasts_1d,
88+
dict(constant_0=[1, 1], smt_0=[1, -1]))
89+
X_1d, contrasts_1d = d_maker(spec_1d, t, level_contrasts=True)
90+
assert_dict_almost_equal(contrasts_1d,
91+
dict(constant_0=1,
92+
smt_0=[1, -1], # main effect
93+
smt_0_0=[1, 0], # level 0, hrf 0
94+
smt_1_0=[0, 1])) # level 1, hrf 0
95+
# Event spec with two factors, one with two levels, another with one
96+
spec_2dc = spec_maker((fac_1, c_fac), ('smt', 'smte'))
97+
X_2dc, contrasts_2dc = d_maker(spec_2dc, t)
98+
assert_almost_equal(X_2dc, np.c_[exp_x_0, exp_x_1])
99+
assert_dict_almost_equal(contrasts_2dc,
100+
{'constant_0': [1, 1],
101+
'smt_0': [1, -1], # main effect
102+
'smt:smte_0': [1, -1], # interaction
103+
})
104+
X_2dc, contrasts_2dc = d_maker(spec_2dc, t, level_contrasts=True)
105+
assert_dict_almost_equal(contrasts_2dc,
106+
{'constant_0': [1, 1],
107+
'smt_0': [1, -1], # main effect
108+
'smt:smte_0': [1, -1], # interaction
109+
'smt_0*smte_1_0': [1, 0], # smt 0, smte 0, hrf 0
110+
'smt_1*smte_1_0': [0, 1], # smt 1, smte 0, hrf 0
111+
})
112+
# Event spec with two factors, both with two levels
113+
spec_2dd = spec_maker((fac_1, fac_2), ('smt', 'smte'))
114+
exp_x_0 = tc_maker((fac_1 == 0) & (fac_2 == 0))
115+
exp_x_1 = tc_maker((fac_1 == 0) & (fac_2 == 1))
116+
exp_x_2 = tc_maker((fac_1 == 1) & (fac_2 == 0))
117+
exp_x_3 = tc_maker((fac_1 == 1) & (fac_2 == 1))
118+
X_2dd, contrasts_2dd = d_maker(spec_2dd, t)
119+
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
120+
exp_cons = {'constant_0': [1, 1, 1, 1],
121+
'smt_0': [1, 1, -1, -1], # main effect fac_1
122+
'smte_0': [1, -1, 1, -1], # main effect fac_2
123+
'smt:smte_0': [1, -1, -1, 1], # interaction
124+
}
125+
assert_dict_almost_equal(contrasts_2dd, exp_cons)
126+
X_2dd, contrasts_2dd = d_maker(spec_2dd, t, level_contrasts=True)
127+
level_cons = exp_cons.copy()
128+
level_cons.update({
129+
'smt_0*smte_0_0': [1, 0, 0, 0], # smt 0, smte 0, hrf 0
130+
'smt_0*smte_1_0': [0, 1, 0, 0], # smt 0, smte 1, hrf 0
131+
'smt_1*smte_0_0': [0, 0, 1, 0], # smt 1, smte 0, hrf 0
132+
'smt_1*smte_1_0': [0, 0, 0, 1], # smt 1, smte 1, hrf 0
133+
})
134+
assert_dict_almost_equal(contrasts_2dd, level_cons)
135+
# Test max order >> 2, no error
136+
X_2dd, contrasts_2dd = d_maker(spec_2dd, t, order=100)
137+
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
138+
assert_dict_almost_equal(contrasts_2dd, exp_cons)
139+
# Test max order = 1
140+
X_2dd, contrasts_2dd = d_maker(spec_2dd, t, order=1)
141+
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
142+
# No interaction
143+
assert_dict_almost_equal(contrasts_2dd,
144+
{'constant_0': [1, 1, 1, 1],
145+
'smt_0': [1, 1, -1, -1], # main effect fac_1
146+
'smte_0': [1, -1, 1, -1], # main effect fac_2
147+
})
148+
# events : test field called "time" is necessary
149+
spec_1d = make_recarray(zip(onsets, fac_1), ('brighteyes', 'smt'))
150+
assert_raises(ValueError, event_design, spec_1d, t)
151+
# blocks : test fields called "start" and "end" are necessary
152+
spec_1d = make_recarray(zip(onsets, offsets, fac_1),
153+
('mister', 'brighteyes', 'smt'))
154+
assert_raises(ValueError, block_design, spec_1d, t)
155+
spec_1d = make_recarray(zip(onsets, offsets, fac_1),
156+
('start', 'brighteyes', 'smt'))
157+
assert_raises(ValueError, block_design, spec_1d, t)
158+
spec_1d = make_recarray(zip(onsets, offsets, fac_1),
159+
('mister', 'end', 'smt'))
160+
assert_raises(ValueError, block_design, spec_1d, t)

0 commit comments

Comments
 (0)