In this notebook we will first implement a simplified version of the DCM model, in order to generate data for subsequent examples.
In [1]:
import numpy
# use a consistent seed so that everyone has the same data
numpy.random.seed(1000)
import os,sys
%matplotlib inline
import matplotlib.pyplot as plt
sys.path.insert(0,'../utils')
from mkdesign import create_design_singlecondition
from graph_utils import show_graph_from_adjmtx
import math
from nipy.modalities.fmri.hemodynamic_models import spm_hrf,compute_regressor
import scipy.interpolate
results_dir = os.path.abspath("../results")
if not os.path.exists(results_dir):
os.mkdir(results_dir)
In [2]:
# first let's build the model without the bilinear influence (aka PPI)
# after http://spm.martinpyka.de/?p=81
nregions=5
z=numpy.zeros(nregions)
# intrinsic connectivity
A=numpy.zeros((z.shape[0],z.shape[0]))
A=numpy.diag(numpy.ones(z.shape[0])*-1)
# add some structure
#A=A + numpy.diag(numpy.ones(z.shape[0]-1),k=-1)
A[2,1]=1
A[3,1]=1
B=numpy.zeros(A.shape)
B[2,0]=1
B[4,0]=1
C=numpy.zeros((z.shape[0],1))
C[0]=1
u=0
print (A)
print (B)
print (C)
# we are assuming a 1 second TR for the resulting data
# but the neural data are at a 1/16 millisecond time resolution
stepsize=.01
tslength=300
timepoints=numpy.arange(0,tslength,stepsize)
# create a blocked design
d,design=create_design_singlecondition(blockiness=1.0,deslength=tslength,blocklength=20,offset=20)
u=scipy.interpolate.griddata(numpy.arange(1,d.shape[0]),d,timepoints,fill_value=0)
def dcm_model(t,z,A,B,C,u):
ut=numpy.abs(timepoints - t).argmin()
return (A.dot(z)+ u[ut]*B.dot(z) + C.dot(u[ut]).T)[0]
def mk_dcm_dataset(timepoints,z,noise_sd):
data=numpy.zeros((len(timepoints),len(z)))
for i in range(1,len(timepoints)):
data[i,:]=data[i-1,:] + dcm_model(timepoints[i],data[i-1,:],A,B,C,u) + numpy.random.randn(len(z))*noise_sd
hrf=spm_hrf(stepsize,oversampling=1)
data_conv=numpy.zeros(data.shape)
for i in range(len(z)):
data_conv[:,i]=numpy.convolve(data[:,i],hrf)[:data.shape[0]]
return data,data_conv
In [3]:
noise_sd=2
data,data_conv=mk_dcm_dataset(timepoints,z,noise_sd)
numpy.savez(os.path.join(results_dir,'dcmdata.npz'),data=data_conv,A=A,B=B,C=C,u=u,d=d,design=design)
In [4]:
plt.subplot(211)
plt.plot(data_conv)
cc=numpy.corrcoef(data_conv.T)
print ('correlation matrix')
print (cc)
from sklearn.covariance import GraphLassoCV
import matplotlib.colors
In [5]:
glasso=GraphLassoCV()
glasso.fit(data_conv)
Out[5]:
In [6]:
from pcor_from_precision import pcor_from_precision
pcor=pcor_from_precision(glasso.precision_)
print ('partial r^2 matrix')
print (pcor**2)
plt.figure(figsize=(10,5))
plt.subplot(141)
plt.imshow(A,interpolation='nearest',norm=matplotlib.colors.Normalize(vmin=-1,vmax=1))
plt.title('A mtx')
plt.subplot(142)
plt.imshow(B,interpolation='nearest',norm=matplotlib.colors.Normalize(vmin=-1,vmax=1))
plt.title('B mtx')
plt.subplot(143)
plt.imshow(cc,interpolation='nearest',norm=matplotlib.colors.Normalize(vmin=-1,vmax=1))
plt.title('correlation')
plt.subplot(144)
plt.imshow(pcor**2,interpolation='nearest',norm=matplotlib.colors.Normalize(vmin=-1,vmax=1))
plt.title('partial correlation')
Out[6]:
In [7]:
gr=show_graph_from_adjmtx(A,B,C)
Now we will run an actual DCM analysis on these data using SPM.
1.Start spm using the command:
spm fmri
2.Click the "Dynamic causal modelling" button on the main panel
3.From the "Action" menu, choose "specify".
The selection menu should automatically go to the directory where the relevant files are (/home/vagrant/fmri-analysis-vm/analysis/connectivity/dcmfiles). Assuming that it does, then first select the SPM.mat file listed in the directory, and click "Done" at the bottom. The SPM.mat file contains a description of the experimental design that DCM needs.
You will be asked for the name: enter "truemodel"
Another file selector will appear, which is asking you to select the ROI files. You should select all 5 of the VOI* files that are listed in the selector. They need to be selected in order; the best way to do this is to hold the Shift key and click the last file in the list, which will select the entire set in order. Click "Done" again.
You will be asked for a number of details about the design, for which you can just accept the default values.
Now we have to specify the three matrices in the SPM model (A,B, and C). You will see a window that asks you to click radio buttons for the active connections. For the A matrix, we need to include both the static and modulated connections from the A and B matrices above. Your entry should look like this:
In [8]:
from IPython.display import Image
Image(filename='DCM_Amtx.png',retina=True)
Out[8]:
Click "done", which will take you to the next window to specify the B and C matrices. Your entries should look like this:
In [9]:
Image(filename='DCM_BCmtx.png',retina=True)
Out[9]:
This will create a file called "DCM_truemodel.mat" in the same directory as the other files.
Now estimate the model by clicking the DCM button, choosing "estimate" from the Action menu, and selecting "DCM_truemodel.mat" from the selector. This will take a few minutes to compute.
Now let's specify the same model, but without any modulated connections. We can do this using MATLAB code rather than doing it by hand.
From the Utils menu, choose "Run M-file". Using the file selector, go one directory up ('..') and select "mk_model_noppi.m". This will specify the model without modulation (i.e. with no nonzero entries in the B matrix) and then estimate it, saving the results to a file called "DCM_noppi.mat'.
Now let's specify the same model, but with two incorrect modulated connections (so that it has the same number of parameters as the true model). Using the "Run M-file" command, select "mk_model_wrongppi.m"
A fundamental part of the DCM workflow is the comparison of models. In this case, we can test for the PPI by comparing the three models that we just created.
Choose "Dynamic causal modelling -> Action -> compare"
This will bring up a batch editor, in which we need to specify several variables.
In [10]:
Image(filename='DCM_BMS.png',retina=True)
Out[10]:
This tells us that the evidence for the first model (the true model) is much greater than that for the other (incorrect) models.
In [ ]: