Skip to content

Commit

Permalink
RF+TST: refactor formula/contrasts to helper func
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@matthew-brett
matthew-brett committed Oct 11, 2015
1 parent fb9bdc1 commit af6cbca
Show file tree
Hide file tree
Showing 2 changed files with 177 additions and 127 deletions.
86 changes: 48 additions & 38 deletions nipy/modalities/fmri/design.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,43 @@ def natural_spline(tvals, knots=None, order=3, intercept=True):
return f.design(tvals, return_float=True)


def _build_formula_contrasts(spec, fields, order):
""" Build formula and contrast in event / block space
Parameters
----------
spec : stuctured array
Structured array containing at least fields listed in `fields`.
fields : sequence of str
Sequence of field names containing names of factors.
order : int
Maximum order of interactions between main effects.
Returns
-------
e_factors : :class:`Formula` instance
Formula for factors given by `fields`
e_contrasts : dict
Dictionary containing contrasts of main effects and interactions
between factors.
"""
e_factors = [Factor(n, np.unique(spec[n])) for n in fields]
e_formula = np.product(e_factors)
e_contrasts = {}
# Add contrasts for factors and factor interactions
max_order = min(len(e_factors), order)
for i in range(1, max_order + 1):
for comb in combinations(zip(fields, e_factors), i):
names = [c[0] for c in comb]
# Collect factors where there is more than one level
fs = [fc.main_effect for fn, fc in comb if len(fc.levels) > 1]
if len(fs) > 0:
e_contrast = np.product(fs).design(spec)
e_contrasts[":".join(names)] = e_contrast
e_contrasts['constant'] = formulae.I.design(spec)
return e_formula, e_contrasts


def event_design(event_spec, t, order=2, hrfs=(glover,),
level_contrasts=False):
""" Create design matrix from event specification `event_spec`
Expand Down Expand Up @@ -131,27 +168,13 @@ def event_design(event_spec, t, order=2, hrfs=(glover,),
itertools.cycle([1])),
('time', '_event_'))
fields = ['_event_']
e_factors = [Factor(n, np.unique(event_spec[n])) for n in fields]
e_formula = np.product(e_factors)

# Design and contrasts in event space
e_formula, e_contrasts = _build_formula_contrasts(
event_spec, fields, order)
# Design and contrasts in block space
# TODO: make it so I don't have to call design twice here
# to get both the contrasts and the e_X matrix as a recarray
e_X = e_formula.design(event_spec)
e_dtype = e_formula.dtype
e_contrasts = {}

# Add contrasts for factors and factor interactions
max_order = min(len(e_factors), order)
for i in range(1, max_order + 1):
for comb in combinations(zip(fields, e_factors), i):
names = [c[0] for c in comb]
# Collect factors where there is more than one level
fs = [fc.main_effect for fn, fc in comb if len(fc.levels) > 1]
if len(fs) > 0:
e_contrast = np.product(fs).design(event_spec)
e_contrasts[":".join(names)] = e_contrast
e_contrasts['constant'] = formulae.I.design(event_spec)

# Now construct the design in time space
t_terms = []
Expand Down Expand Up @@ -229,40 +252,27 @@ def block_design(block_spec, t, order=2, hrfs=(glover,),
raise ValueError('expecting fields called "start" and "end"')
fields.pop(fields.index('start'))
fields.pop(fields.index('end'))
e_factors = [Factor(n, np.unique(block_spec[n])) for n in fields]
e_formula = np.product(e_factors)
e_contrasts = {}

if len(e_factors) > 1:
for i in range(1, order+1):
for comb in combinations(zip(fields, e_factors), i):
names = [c[0] for c in comb]
fs = [c[1].main_effect for c in comb]
e_contrasts[":".join(names)] = np.product(fs).design(block_spec)
else:
# only one factor, produce the main effect
field = fields[0]
factor = e_factors[0]
e_contrasts[field + '_effect'] = factor.main_effect.design(event_spec)

e_contrasts['constant'] = formulae.I.design(block_spec)

if len(fields) == 0: # No factors specified, make generic block
block_spec = make_recarray(zip(block_spec['start'],
block_spec['end'],
itertools.cycle([1])),
('start', 'end', '_block_'))
fields = ['_block_']
e_formula, e_contrasts = _build_formula_contrasts(
block_spec, fields, order)
# Design and contrasts in block space
# TODO: make it so I don't have to call design twice here
# to get both the contrasts and the e_X matrix as a recarray

e_X = e_formula.design(block_spec)
e_dtype = e_formula.dtype

# Now construct the design in time space

block_times = np.array([(s,e) for s, e in zip(block_spec['start'],
block_spec['end'])])
convolution_interval = (block_times.min() - convolution_padding,
block_times.max() + convolution_padding)

t_terms = []
t_names = []
t_contrasts = {}
for l, h in enumerate(hrfs):
for n in e_dtype.names:
Expand Down
218 changes: 129 additions & 89 deletions nipy/modalities/fmri/tests/test_design.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@
import numpy as np

from ..design import event_design, block_design
from ..utils import events, lambdify_t
from ..utils import (events, lambdify_t, T, convolve_functions,
blocks)
from ..hrf import glover
from nipy.algorithms.statistics.formula import make_recarray

Expand All @@ -26,95 +27,134 @@ def test_event_design():
# Test event design helper function
# An event design with one event type
onsets = np.array([0, 20, 40, 60])
durations = np.array([2, 3, 4, 5])
offsets = onsets + durations
c_fac = np.array([1, 1, 1, 1]) # constant factor (one level)
fac_1 = np.array([0, 1, 0, 1]) # factor 1, two levels
fac_2 = np.array([0, 0, 1, 1]) # factor 2, two levels
t = np.arange(0, 100, 2)
# Event spec with no event factor -> single column design, no contrasts
event_spec_0 = make_recarray(onsets, ('time',))
X_0, contrasts_0 = event_design(event_spec_0, t)
exp_x_0 = lambdify_t(events(onsets, f=glover))(t)
assert_almost_equal(X_0, exp_x_0)
assert_dict_almost_equal(contrasts_0, dict(constant_0=1))
X_0, contrasts_0 = event_design(event_spec_0, t, level_contrasts=True)
assert_almost_equal(X_0, exp_x_0)
assert_dict_almost_equal(contrasts_0, dict(constant_0=1, _event__1_0=1))
# Event spec with single factor, but only one level
event_spec_1c = make_recarray(zip(onsets, c_fac), ('time', 'smt'))
X_1c, contrasts_1c = event_design(event_spec_1c, t)
assert_almost_equal(X_1c, exp_x_0)
assert_dict_almost_equal(contrasts_1c, dict(constant_0=1))
X_1c, contrasts_1c = event_design(event_spec_1c, t, level_contrasts=True)
assert_dict_almost_equal(contrasts_1c, dict(constant_0=1, smt_1_0=1))
# Event spec with single factor, two levels
event_spec_1d = make_recarray(zip(onsets, fac_1), ('time', 'smt'))
exp_x_0 = lambdify_t(events(onsets[fac_1 == 0], f=glover))(t)
exp_x_1 = lambdify_t(events(onsets[fac_1 == 1], f=glover))(t)
X_1d, contrasts_1d = event_design(event_spec_1d, t)
assert_almost_equal(X_1d, np.c_[exp_x_0, exp_x_1])
assert_dict_almost_equal(contrasts_1d,
dict(constant_0=[1, 1], smt_0=[1, -1]))
X_1d, contrasts_1d = event_design(event_spec_1d, t, level_contrasts=True)
assert_dict_almost_equal(contrasts_1d,
dict(constant_0=1,
smt_0=[1, -1], # main effect
smt_0_0=[1, 0], # level 0, hrf 0
smt_1_0=[0, 1])) # level 1, hrf 0
# Event spec with two factors, one with two levels, another with one
event_spec_2dc = make_recarray(zip(onsets, fac_1, c_fac),
('time', 'smt', 'smte'))
X_2dc, contrasts_2dc = event_design(event_spec_2dc, t)
assert_almost_equal(X_2dc, np.c_[exp_x_0, exp_x_1])
assert_dict_almost_equal(contrasts_2dc,
{'constant_0': [1, 1],
'smt_0': [1, -1], # main effect
'smt:smte_0': [1, -1], # interaction
})
X_2dc, contrasts_2dc = event_design(event_spec_2dc, t, level_contrasts=True)
assert_dict_almost_equal(contrasts_2dc,
{'constant_0': [1, 1],
'smt_0': [1, -1], # main effect
'smt:smte_0': [1, -1], # interaction
'smt_0*smte_1_0': [1, 0], # smt 0, smte 0, hrf 0
'smt_1*smte_1_0': [0, 1], # smt 1, smte 0, hrf 0
})
# Event spec with two factors, both with two levels
event_spec_2dd = make_recarray(zip(onsets, fac_1, fac_2),
('time', 'smt', 'smte'))
exp_x_0 = lambdify_t(events(onsets[(fac_1 == 0) & (fac_2 == 0)], f=glover))(t)
exp_x_1 = lambdify_t(events(onsets[(fac_1 == 0) & (fac_2 == 1)], f=glover))(t)
exp_x_2 = lambdify_t(events(onsets[(fac_1 == 1) & (fac_2 == 0)], f=glover))(t)
exp_x_3 = lambdify_t(events(onsets[(fac_1 == 1) & (fac_2 == 1)], f=glover))(t)
X_2dd, contrasts_2dd = event_design(event_spec_2dd, t)
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
exp_cons = {'constant_0': [1, 1, 1, 1],
'smt_0': [1, 1, -1, -1], # main effect fac_1
'smte_0': [1, -1, 1, -1], # main effect fac_2
'smt:smte_0': [1, -1, -1, 1], # interaction
}
assert_dict_almost_equal(contrasts_2dd, exp_cons)
X_2dd, contrasts_2dd = event_design(event_spec_2dd, t, level_contrasts=True)
level_cons = exp_cons.copy()
level_cons.update({
'smt_0*smte_0_0': [1, 0, 0, 0], # smt 0, smte 0, hrf 0
'smt_0*smte_1_0': [0, 1, 0, 0], # smt 0, smte 1, hrf 0
'smt_1*smte_0_0': [0, 0, 1, 0], # smt 1, smte 0, hrf 0
'smt_1*smte_1_0': [0, 0, 0, 1], # smt 1, smte 1, hrf 0
})
assert_dict_almost_equal(contrasts_2dd, level_cons)
# Test max order >> 2, no error
X_2dd, contrasts_2dd = event_design(event_spec_2dd, t, order=100)
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
assert_dict_almost_equal(contrasts_2dd, exp_cons)
# Test max order = 1
X_2dd, contrasts_2dd = event_design(event_spec_2dd, t, order=1)
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
# No interaction
assert_dict_almost_equal(contrasts_2dd,
{'constant_0': [1, 1, 1, 1],
'smt_0': [1, 1, -1, -1], # main effect fac_1
'smte_0': [1, -1, 1, -1], # main effect fac_2
})
# Test field called "time" is necessary
event_spec_1d = make_recarray(zip(onsets, fac_1), ('brighteyes', 'smt'))
assert_raises(ValueError, event_design, event_spec_1d, t)

def mk_ev_spec(factors, names):
names = ('time',) + tuple(names)
if len(factors) == 0:
return make_recarray(onsets, names)
return make_recarray(zip(onsets, *factors), names)

def mk_blk_spec(factors, names):
names = ('start', 'end') + tuple(names)
return make_recarray(zip(onsets, offsets, *factors), names)

def mk_ev_tc(ev_inds):
# Make real time course for given event onsets
return lambdify_t(events(onsets[ev_inds], f=glover))(t)

def mk_blk_tc(ev_inds):
# Make real time course for block onset / offsets
B = blocks(zip(onsets[ev_inds], offsets[ev_inds]))
term = convolve_functions(B, glover(T),
(-5, 70), # step func support
(0, 30.), # conv kernel support
0.02) # dt
return lambdify_t(term)(t)

for d_maker, spec_maker, tc_maker, null_name in (
(event_design, mk_ev_spec, mk_ev_tc, '_event_'),
(block_design, mk_blk_spec, mk_blk_tc, '_block_')):
# Event spec with no event factor -> single column design, no contrasts
spec_0 = spec_maker((), ())
X_0, contrasts_0 = d_maker(spec_0, t)
exp_x_0 = tc_maker(onsets==onsets)
assert_almost_equal(X_0, exp_x_0)
assert_dict_almost_equal(contrasts_0, dict(constant_0=1))
X_0, contrasts_0 = d_maker(spec_0, t, level_contrasts=True)
assert_almost_equal(X_0, exp_x_0)
assert_dict_almost_equal(contrasts_0,
{'constant_0': 1,
null_name + '_1_0': 1})
# Event spec with single factor, but only one level
spec_1c = spec_maker((c_fac,), ('smt',))
X_1c, contrasts_1c = d_maker(spec_1c, t)
assert_almost_equal(X_1c, exp_x_0)
assert_dict_almost_equal(contrasts_1c, dict(constant_0=1))
X_1c, contrasts_1c = d_maker(spec_1c, t, level_contrasts=True)
assert_dict_almost_equal(contrasts_1c, dict(constant_0=1, smt_1_0=1))
# Event spec with single factor, two levels
spec_1d = spec_maker((fac_1,), ('smt',))
exp_x_0 = tc_maker(fac_1 == 0)
exp_x_1 = tc_maker(fac_1 == 1)
X_1d, contrasts_1d = d_maker(spec_1d, t)
assert_almost_equal(X_1d, np.c_[exp_x_0, exp_x_1])
assert_dict_almost_equal(contrasts_1d,
dict(constant_0=[1, 1], smt_0=[1, -1]))
X_1d, contrasts_1d = d_maker(spec_1d, t, level_contrasts=True)
assert_dict_almost_equal(contrasts_1d,
dict(constant_0=1,
smt_0=[1, -1], # main effect
smt_0_0=[1, 0], # level 0, hrf 0
smt_1_0=[0, 1])) # level 1, hrf 0
# Event spec with two factors, one with two levels, another with one
spec_2dc = spec_maker((fac_1, c_fac), ('smt', 'smte'))
X_2dc, contrasts_2dc = d_maker(spec_2dc, t)
assert_almost_equal(X_2dc, np.c_[exp_x_0, exp_x_1])
assert_dict_almost_equal(contrasts_2dc,
{'constant_0': [1, 1],
'smt_0': [1, -1], # main effect
'smt:smte_0': [1, -1], # interaction
})
X_2dc, contrasts_2dc = d_maker(spec_2dc, t, level_contrasts=True)
assert_dict_almost_equal(contrasts_2dc,
{'constant_0': [1, 1],
'smt_0': [1, -1], # main effect
'smt:smte_0': [1, -1], # interaction
'smt_0*smte_1_0': [1, 0], # smt 0, smte 0, hrf 0
'smt_1*smte_1_0': [0, 1], # smt 1, smte 0, hrf 0
})
# Event spec with two factors, both with two levels
spec_2dd = spec_maker((fac_1, fac_2), ('smt', 'smte'))
exp_x_0 = tc_maker((fac_1 == 0) & (fac_2 == 0))
exp_x_1 = tc_maker((fac_1 == 0) & (fac_2 == 1))
exp_x_2 = tc_maker((fac_1 == 1) & (fac_2 == 0))
exp_x_3 = tc_maker((fac_1 == 1) & (fac_2 == 1))
X_2dd, contrasts_2dd = d_maker(spec_2dd, t)
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
exp_cons = {'constant_0': [1, 1, 1, 1],
'smt_0': [1, 1, -1, -1], # main effect fac_1
'smte_0': [1, -1, 1, -1], # main effect fac_2
'smt:smte_0': [1, -1, -1, 1], # interaction
}
assert_dict_almost_equal(contrasts_2dd, exp_cons)
X_2dd, contrasts_2dd = d_maker(spec_2dd, t, level_contrasts=True)
level_cons = exp_cons.copy()
level_cons.update({
'smt_0*smte_0_0': [1, 0, 0, 0], # smt 0, smte 0, hrf 0
'smt_0*smte_1_0': [0, 1, 0, 0], # smt 0, smte 1, hrf 0
'smt_1*smte_0_0': [0, 0, 1, 0], # smt 1, smte 0, hrf 0
'smt_1*smte_1_0': [0, 0, 0, 1], # smt 1, smte 1, hrf 0
})
assert_dict_almost_equal(contrasts_2dd, level_cons)
# Test max order >> 2, no error
X_2dd, contrasts_2dd = d_maker(spec_2dd, t, order=100)
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
assert_dict_almost_equal(contrasts_2dd, exp_cons)
# Test max order = 1
X_2dd, contrasts_2dd = d_maker(spec_2dd, t, order=1)
assert_almost_equal(X_2dd, np.c_[exp_x_0, exp_x_1, exp_x_2, exp_x_3])
# No interaction
assert_dict_almost_equal(contrasts_2dd,
{'constant_0': [1, 1, 1, 1],
'smt_0': [1, 1, -1, -1], # main effect fac_1
'smte_0': [1, -1, 1, -1], # main effect fac_2
})
# events : test field called "time" is necessary
spec_1d = make_recarray(zip(onsets, fac_1), ('brighteyes', 'smt'))
assert_raises(ValueError, event_design, spec_1d, t)
# blocks : test fields called "start" and "end" are necessary
spec_1d = make_recarray(zip(onsets, offsets, fac_1),
('mister', 'brighteyes', 'smt'))
assert_raises(ValueError, block_design, spec_1d, t)
spec_1d = make_recarray(zip(onsets, offsets, fac_1),
('start', 'brighteyes', 'smt'))
assert_raises(ValueError, block_design, spec_1d, t)
spec_1d = make_recarray(zip(onsets, offsets, fac_1),
('mister', 'end', 'smt'))
assert_raises(ValueError, block_design, spec_1d, t)

0 comments on commit af6cbca

Please sign in to comment.