How to specify a DCM/rDCM

Simple DCM

In this toy example we show how to create a simple linear DCM with 3 regions. First we define the network architecture (presence and absence of connections).

using RegressionDynamicCausalModeling
using Random

a = BitMatrix([1 0 1; 1 1 0; 1 0 1])
c = BitMatrix([0 0; 1 0; 0 1])

3×2 BitMatrix:
 0  0
 1  0
 0  1

Number of datapoints per region

scans = 300

300

Specify the information about driving inputs

off = zeros(100)
on = ones(100)
u = [off; on; off; on; off]
u = [u u]

u_dt = 1 / 32 # sampling interval
name = ["input1", "input2"]

U = InputU(u, u_dt, name)

U (Input)
   u:  500x2 matrix
   dt: 0.03125s
   names: input1,...,input2

Next, we specify information about the BOLD signal

y = nothing
y_dt = 0.5 # sampling interval
name = ["region1", "region2", "region3"]

Y = BoldY(y, y_dt, name)

Y (BOLD signal)
   y: empty
   dt: 0.5s
   names: empty

The last specification we need are the actual parameter values of the A and C matrix. In this example we generate values by sampling from the prior.

Ep = rDCM.TrueParamLinear(a, c; sample=true, rng=MersenneTwister(42))

Ep (True parameters)
    A: 3 x 3 matrix
    C: 3 x 2 matrix
    transit: 0.00 ... 0.00
    decay:   0.00 ... 0.00
    epsilon: 0.00

We can now construct our linear DCM and use it for subsequent analysis

dcm = LinearDCM(a, c, scans, size(a, 1), U, Y, Ep)

Linear DCM
a:     3x3 matrix
c:     3x2 matrix
scans: 300
nr:    3
---------------------------------------------
U (Input)
   u:  500x2 matrix
   dt: 0.03125s
   names: input1,...,input2
---------------------------------------------
Y (BOLD signal)
   y: empty
   dt: 0.5s
   names: empty
---------------------------------------------
Ep (True parameters)
    A: 3 x 3 matrix
    C: 3 x 2 matrix
    transit: 0.00 ... 0.00
    decay:   0.00 ... 0.00
    epsilon: 0.00

Confounds

It is possible to specify confounds. The following example shows how this can be done.

X0 = zeros(500, 3) # specify confound matrix (fist dimension is time, second dimension is number of confounds)
conf_names = ["conf1", "conf2", "conf3"]
conf = Confound(X0, conf_names)

Confounds
   X0:    500x3 matrix
   names: conf1,...,conf3

Construct DCM and convert to rDCM structure

dcm = LinearDCM(a, c, scans, size(a, 1), U, Y, Ep, conf)
rdcm = RigidRdcm(dcm)

rigid rDCM
a:     3x3 matrix
c:     3x2 matrix
scans: 300
nr:    3
HRF:   500 element vector
---------------------------------------------
U (Input)
   u:  500x2 matrix
   dt: 0.03125s
   names: input1,...,input2
---------------------------------------------
Y (BOLD signal)
   y: empty
   dt: 0.5s
   names: empty
---------------------------------------------
Ep (True parameters)
    A: 3 x 3 matrix
    C: 3 x 2 matrix
    transit: 0.00 ... 0.00
    decay:   0.00 ... 0.00
    epsilon: 0.00
---------------------------------------------
Confounds
   X0:    500x3 matrix
   names: conf1,...,conf3

This page was generated using Literate.jl.