<?xml version="1.0" encoding="utf-8"?><!DOCTYPE article  PUBLIC '-//OASIS//DTD DocBook XML V4.4//EN'  'http://www.docbook.org/xml/4.4/docbookx.dtd'><article><articleinfo><title>CodeSnippet</title><revhistory><revision><revnumber>2</revnumber><date>2013-03-08 10:28:25</date><authorinitials>localhost</authorinitials><revremark>converted to 1.6 markup</revremark></revision><revision><revnumber>1</revnumber><date>2007-07-27 17:11:31</date><authorinitials>devel06.mrc-cbu.cam.ac.uk</authorinitials></revision></revhistory></articleinfo><screen><![CDATA[''' Run models on ROI data '''
]]><![CDATA[
import os
]]><![CDATA[
import numpy as N
from scipy import io as sio
import scipy.sandbox.models as SSM
import scipy.interpolate as SI
]]><![CDATA[
from neuroimaging.modalities.fmri import protocol
from neuroimaging.modalities.fmri import hrf
from neuroimaging.modalities.fmri import functions
]]><![CDATA[
data_file = os.path.join('..', 'data', 'roi_data.mat')
roi_data = sio.loadmat(data_file)['roi_data']
]]><![CDATA[
# Specify type of linear model to estimate
model_type = SSM.regression.ar_model
ar_rho = 3
ar_niter = 6
]]><![CDATA[
# Specify drift term for all models
drift_window = (0,256)
drift_df = 7
drift = functions.SplineConfound(df=drift_df,
                                 window=drift_window)
drift_term = protocol.ExperimentalQuantitative(
    'spline_drift',
    drift)
]]><![CDATA[
analyses = {}
for sb in range(roi_data.shape[0]):
    analyses[sb] = {}
    for ss in range(roi_data.shape[1]):
        data = roi_data[sb, ss]
        onsets = data.ons
        tr = data.TR
        y = data.y
        movements = data.movements
        # Adust onsets to seconds
        onsets *= tr
        # Create time vector
        T = N.arange(len(y)) * tr
        # Make event factor
        flashes = zip(['flash']*len(onsets), onsets)
        flash_factor = protocol.ExperimentalFactor('flash',
                                                   flashes,
                                                   delta=True)
        flash_factor.convolve(hrf.canonical)
        # Add movements by providing interpolator
        moves_interpolator = functions.InterpolatedConfound(T, movements.T)
        movement_term = protocol.ExperimentalQuantitative(
            'movement',
            moves_interpolator)
        # Add factor, moves and drift term to model
        formula = flash_factor + movement_term + drift_term
        # Calculate design from model and time vector
        design = formula(time=T).T
        # Estimate design on time series data
        model = model_type(design, rho=ar_rho)
        result = model.iterative_fit(y, niter=ar_niter)
        analyses[sb][ss] = {'model': model,
                            'result': result}]]></screen></article>