In this notebook, we will work through a simulation of psychophysiological interaction



In [ ]:

    
import os,sys
import numpy
%matplotlib inline
import matplotlib.pyplot as plt
sys.path.insert(0,'../')
from utils.mkdesign import create_design_singlecondition
from nipy.modalities.fmri.hemodynamic_models import spm_hrf,compute_regressor
from utils.make_data import make_continuous_data
from utils.corr_timeseries.exact_switching import generate_data as generate_exact
from statsmodels.tsa.arima_process import arma_generate_sample
import scipy.stats

Load the data generated using the DCM forward model. In this model, there should be a significant PPI between roi 0 and rois 2 and 4 (see the B matrix in the DCM notebook)



In [ ]:

    
data_conv=numpy.load('dcmdata.npy')
# downsample to 1 second TR
data=data_conv[range(0,data_conv.shape[0],100)]
ntp=data.shape[0]

# create a blocked design
d,design=create_design_singlecondition(blockiness=1.0,deslength=ntp,blocklength=20,offset=20)

regressor,_=compute_regressor(design,'spm',numpy.arange(0,ntp))



In [ ]:

    
plt.plot(data)

Set up the PPI model, using ROI 0 as the seed



In [ ]:

    
seed=0
X=numpy.vstack((regressor[:,0],data[:,seed],regressor[:,0]*data[:,seed],numpy.ones(data.shape[0]))).T
hat_mtx=numpy.linalg.inv(X.T.dot(X)).dot(X.T)

for i in range(data.shape[1]):
    if i==seed:
        continue
    beta_hat=hat_mtx.dot(data[:,i])
    resid=data[:,i] - X.dot(beta_hat)
    sigma2hat=(resid.dot(resid))/(X.shape[0] - X.shape[1])
    c=numpy.array([0,0,1,0])  # contrast for PPI
    t=c.dot(beta_hat)/numpy.sqrt(c.dot(numpy.linalg.inv(X.T.dot(X)).dot(c))*sigma2hat)
    print 'ROI %d:'%i, t, 1.0 - scipy.stats.t.cdf(t,X.shape[0] - X.shape[1])

Let's plot the relation between the ROIs as a function of the task



In [ ]:

    
on_tp=numpy.where(regressor>0.9)[0]
off_tp=numpy.where(regressor<0.01)[0]

plt.scatter(data[on_tp,0],data[on_tp,2])
fit = np.polyfit(data[on_tp,0],data[on_tp,2],1)
plt.plot()
plt.scatter(data[off_tp,0],data[off_tp,2],color='red')



In [ ]:

    
foo=regressor[:,0]*data[:,0]



In [ ]:

    
foo.shape



In [ ]:

    
data[:,0].shape



In [ ]: