DGM-Simulations

Packages and main variables

Install required packages

# install.packages("rmarkdown")
# install.packages("multdyn")
# install.packages("R.matlab")
# install.packages("cowplot")
# install.packages("png")
# install.packages("testit")

Load libraries

library(DGM)
library(R.matlab)
library(testit)
library(ggplot2)
library(cowplot)
library(reshape2)
library(png)
library(grid)

Main variables

N=50 # Number of simulated subjects/datasets
Nn=5 # Number of nodes
PATH_HOME = "/home/simon"
PATH = file.path(PATH_HOME, "Dropbox", "Data", "DGM-Sim")  # Project path
PATH_FIG  = file.path(PATH, 'figures') # path where figures will be stored
PATH_RES  = file.path(PATH, 'results') # path where results will be stored
PATH_TS = file.path(PATH, 'data', 'sim', 'timeseries') # path where time series data is
PATH_NET = file.path(PATH, 'data', 'sim', 'nets') # path where network data is
Sys.setenv(R_PATH_TS = PATH_TS)

Get Sim1 and Sim22 from FMRIB

if [ -f ${R_PATH_TS}/sim1.mat ]; then
  echo Found sim1 and sim22.
else
  echo Downloading sim1 and sim22...
  wget http://www.fmrib.ox.ac.uk/datasets/netsim/sims.tar.gz -P ${R_PATH_TS} >/dev/null 2>&1
  tar zxvf ${R_PATH_TS}/sims.tar.gz -C ${R_PATH_TS} sim1.mat sim22.mat
  rm ${R_PATH_TS}/sims.tar.gz
fi

Found sim1 and sim22.

Loading time series data

# Downloaded from http://www.fmrib.ox.ac.uk/datasets/netsim/
S=200 # No of samples for Sim1 and Sim22
d = readMat(file.path(PATH_TS,'sim1.mat'))
ts.sim1 = reshapeTs(d$ts,N,S)
d = readMat(file.path(PATH_TS,'sim22.mat'))
ts.sim22 = reshapeTs(d$ts,N,S)
ts.int0 = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF1_Mod1_Inj0_F1.mat'))$gfy2s
ts.int1 = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F13.mat'))$gfy2s
ts.int2 = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F16.mat'))$gfy2s
ts.int3 = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F20.mat'))$gfy2s
ts.int4 = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F24.mat'))$gfy2s
ts.int5 = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F28.mat'))$gfy2s
ts.int6 = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F32.mat'))$gfy2s
# Very long 60 min simulation
ts.long = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF4_Mod1_Inj0_F1_60min.mat'))$gfy2s
# noise
ts.noise = readMat(file.path(PATH_TS,'Nn5_TR2_Noise1_HRF4_Mod1_Inj0_F1_noise.mat'))$gfy2s

Plot timeseries of random subject

t = 1:50 # interval to plot
set.seed(1980)
s = sample(N,1) # random subject
vn = c("time", "node")
d=list()
d[[1]] = melt(ts.sim1[t,,s], varnames = vn)
d[[2]] = melt(ts.sim22[t,,s],varnames = vn)
d[[3]] = melt(ts.int0[t,,s], varnames = vn)
d[[4]] = melt(ts.int1[t,,s], varnames = vn)
d[[5]] = melt(ts.int2[t,,s], varnames = vn)
d[[6]] = melt(ts.int3[t,,s], varnames = vn)
d[[7]] = melt(ts.int4[t,,s], varnames = vn)
d[[8]] = melt(ts.int5[t,,s], varnames = vn)
d[[9]] = melt(ts.int6[t,,s], varnames = vn)
d[[10]] = melt(ts.long[t,,s], varnames = vn)
d[[11]] = melt(ts.noise[t,,s], varnames = vn)
p=list()
str_int = c("sim1", "sim22", "int0", "int1", "int2", "int3", "int4", "int5", "int6" , "long", "noise")
for (i in 1:length(d)) {
  p[[i]] = ggplot(d[[i]], aes(x = time, y = value, group=node, color=as.factor(node))) + geom_line() +
    theme_minimal() + ggtitle(str_int[i]) + scale_color_discrete(name = "node")
}
plot_grid(plotlist = p, ncol = 2, nrow = 6, rel_widths = c(1, 1))

Load time series data (single spike)

ts.ssint0 = readMat(file.path(PATH_TS, "SingleSpike_Nn5_TR01_Noise01_HRF1_Mod1_Inj0_F1.mat"))$gytrue
ts.ssint1 = readMat(file.path(PATH_TS, "SingleSpike_Nn5_TR01_Noise01_HRF1_Mod1_Inj1_F13.mat"))$gytrue
ts.ssint2 = readMat(file.path(PATH_TS, "SingleSpike_Nn5_TR01_Noise01_HRF1_Mod1_Inj1_F16.mat"))$gytrue
ts.ssint3 = readMat(file.path(PATH_TS, "SingleSpike_Nn5_TR01_Noise01_HRF1_Mod1_Inj1_F20.mat"))$gytrue
ts.ssint4 = readMat(file.path(PATH_TS, "SingleSpike_Nn5_TR01_Noise01_HRF1_Mod1_Inj1_F24.mat"))$gytrue
ts.ssint5 = readMat(file.path(PATH_TS, "SingleSpike_Nn5_TR01_Noise01_HRF1_Mod1_Inj1_F28.mat"))$gytrue
ts.ssint6 = readMat(file.path(PATH_TS, "SingleSpike_Nn5_TR01_Noise01_HRF1_Mod1_Inj1_F32.mat"))$gytrue

Estimate networks (example for a single simulation data set)

# for (s in 1:N) {
#  s=subject(scaleTs(ts.sim1[,,s]), id=sprintf("Id_%03d", s), 
#            path = file.path(PATH_NET, "sim1"))
# }
# 
# for (s in 1:N) {
#   s=subject(scaleTs(ts.noise[,,s]), id=sprintf("Id_%03d", s),
#             path = file.path(PATH_NET, "Nn5_TR2_Noise1_HRF4_Mod1_Inj0_F1_noise"))
# }
# 
# for (s in 1:N) {
#   s=subject(scaleTs(ts.noise[,,s]), id=sprintf("Id_%03d", s),
#             path = file.path(PATH_NET, "Nn5_TR3_Noise01_HRF4_Mod1_Inj0_F1"))
# }

Generate true network

atrue=array(0,dim=c(5,5))
atrue[1,2] = atrue[2,3] = atrue[3,4] = atrue[4,5] = atrue[1,5] = 1
btrue = atrue==1
example=atrue
example[2,1]=1
example[4,5]=0
example[5,4]=1
p1=gplotMat(atrue, title = "true network", hasColMap = F)
p2=gplotMat(t(atrue), title = "inverse directionality", hasColMap = F)
p3=gplotMat(t(atrue)+atrue, title = "bidirectional", hasColMap = F)
p4=gplotMat(example, title = "example", hasColMap = F)
plot_grid(p1, p2, p3, p4, ncol = 4, nrow = 1)

ggsave(path = PATH_FIG, "trueNetAndVariants.png")

Saving 7.2 x 2 in image

#perf(atrue, atrue)
#perf(t(atrue), atrue)
#perf(t(atrue)+atrue, atrue)
perf(example, atrue)

$subj
     tpr       spc       ppv       npv       fpr fnr       fdr  acc
[1,] 0.8 0.8666667 0.6666667 0.9285714 0.1333333 0.2 0.3333333 0.85

$cases
     TP FP FN TN
[1,]  4  2  1 13

$tpr
[1] 0.8

$spc
[1] 0.8666667

$acc
[1] 0.85

$ppv
[1] 0.6666667

Computation benchmarks

Commented code was run on a execution node Intel Xeon CPU E5-2630 v2 @ 2.60GHz with R 3.4.0

# n=13
# t=1200
# k=3:n
# 
# time=rep(NA,1,length(k))
# X=array(rnorm(t*n), dim=c(t,n))
# 
# c=1;
# for (i in k) {
#   time[c]=system.time(exhaustive.search(X[,1:i],1))[3]
#   c=c+1
# }
# # Quick bench 8 nodes
# X=array(rnorm(1200*8), dim=c(1200,8))
# system.time(exhaustive.search(X,1))[3]
# execution time values from buster
k=3:13
time = c(0.569, 1.117, 2.357, 5.081, 11.103, 24.035, 51.979, 112.249,
         240.239, 505.230, 1098.665)
time = c(0.257, 0.518, 0.920, 2.013, 4.243, 9.263, 19.855, 42.651, 91.377,
         194.460, 399.468)
fit = lm(log(time) ~ k)
# plot(k, time, pch=16)
j=c(15,20,25)
r=exp(fit$coefficients[1] + fit$coefficients[2]*j)
nodes=c(k,j)
time=c(time,r)
x=rbind(nodes, time)
# print(x, digits = 2)
f=c(rep(1,8), rep(60,4),60^2, 60^2*24)
x[2,]=x[2,]/f
print(x, digits = 2)

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
nodes 3.00 4.00 5.00    6  7.0  8.0    9   10 11.0  12.0  13.0    15    20    25
time  0.26 0.52 0.92    2  4.2  9.3   20   43  1.5   3.2   6.7    29    20    35

Quick Bench of 8-node networks

X=array(rnorm(1200*8), dim=c(1200,8))
# system.time(exhaustive.search(X,1))[3]

Figure 1: True network and correlation matrix

comput=data.frame(nodes = as.factor(nodes), time= time)
p5 = ggplot(data=comput, aes(x=nodes, y=time^(1/3))) + 
  geom_bar(stat="identity", fill="steelblue") +
  geom_text(aes(label=c("0.3\nsec","0.5\nsec","0.9\nsec","2.0\nsec","4.2\nsec",
                        "9.3\nsec","20\nsec", "43\nsec","1.5\nmin","3.2\nmin","6.7\nmin",
                        "29\nmin","20\nhrs","35\ndays")), size=2.5, vjust=-0.3) +
  theme_minimal() + ylab(expression('time s'^(1/3))) + ylim(c(0, 210)) + xlab("network size")
img = readPNG(file.path(PATH_FIG, "fig-truenet-page001.png"))
g = rasterGrob(img, interpolate=T)
p1 = ggplot() + annotation_custom(g) + ggtitle('Simulated\n5-node network')
p2 = gplotMat(rmna(btrue), title='5-node\nnetwork', hasColMap = F)
p3 = gplotMat(rmdiag(corTs(ts.sim22)), title='Dynamic', barWidth = 0.2,
              colMapLabel = expression("Pearson\'s"~italic(r)), lim = c(0, 0.5)) + xlab("Node") + ylab("Node")
p4 = gplotMat(rmdiag(corTs(ts.sim1)), title='Stationary', barWidth = 0.2,
              colMapLabel = expression("Pearson\'s"~italic(r)), lim = c(0, 0.5)) +  xlab("Node") + ylab("Node")
a = plot_grid(p1, p2, ncol=2, nrow = 1, rel_widths = c(1, 0.8), labels="A")
c = plot_grid(p5, ncol=1, labels = "C")
left = plot_grid(a, c, ncol=1,  rel_heights = c(0.9, 1))
right = plot_grid(p3, p4, ncol=1, nrow=2, labels = "B")
plot_grid(left, right, ncol=2, rel_widths = c(1, 0.85))

ggsave(path = PATH_FIG, "Fig1.png")

Signal standard deviation

SD_sim22 = SD_int0 = SD_sim1 =array(NA, dim=c(N,Nn))
for (i in 1:N) {
  SD_sim22[i,]= apply(ts.sim22[,,i], 2, sd)
  SD_int0[i,] = apply(ts.int0[,,i], 2, sd)
  SD_sim1[i,] = apply(ts.sim1[,,i], 2, sd)
}
x=t(array(c(colMeans(SD_sim22), colMeans(SD_int0), colMeans(SD_sim1)), dim=c(5,3)))
colnames(x)=c("node1", "node2", "node3", "node4", "node5")
rownames(x)=c("Sim22", "int0", "Sim1")
print(x)

         node1    node2    node3    node4    node5
Sim22 2.026333 1.419497 1.196420 1.088113 1.557278
int0  2.002263 1.524760 1.248987 1.137556 1.632403
Sim1  2.203949 2.281024 2.271040 2.263299 2.336683

Signal SD (mean across subjects). Variability decreases from node 1 to node 4 with node 5 having higher variability. Consistant with simulation 22.

Global mean of SD

rowMeans(x)

   Sim22     int0     Sim1 
1.457528 1.509194 2.271199 

Pearson’s correlations of the nodes

R=array(NA,dim=c(10,Nn)) # Sim. dataset x nodes
R[1,] = corTs(ts.sim1)[btrue]
R[2,] = corTs(ts.sim22)[btrue]
R[3,]= corTs(ts.long)[btrue]
R[4,] = corTs(ts.int0)[btrue]
R[5,] = corTs(ts.int1)[btrue]
R[6,] = corTs(ts.int2)[btrue]
R[7,] = corTs(ts.int3)[btrue]
R[8,] = corTs(ts.int4)[btrue]
R[9,] = corTs(ts.int5)[btrue]
R[10,] = corTs(ts.int6)[btrue]
idx = c(1,2,3,5,4) # move connection 1->5 to last column
colnames(R)=c("1->2", "2->3", "3->4", "4->5", "1->5")
rownames(R)=c("sim1", "sim22", "long", "int0", "int1", "int2",
              "int3", "int4", "int5", "int6")
print(R[,idx])

           1->2      2->3      3->4      1->5      4->5
sim1  0.3054908 0.3488552 0.3214245 0.3319976 0.3011353
sim22 0.3978376 0.3936898 0.2927097 0.2481423 0.3381125
long  0.3876135 0.3294120 0.3242838 0.2102027 0.3633495
int0  0.4211651 0.3543212 0.3702699 0.2317834 0.3857967
int1  0.4006303 0.3575501 0.3671319 0.2227266 0.3775324
int2  0.3889032 0.3558315 0.3645721 0.2147245 0.3614009
int3  0.3684931 0.3532474 0.3571739 0.2006717 0.3367386
int4  0.3466638 0.3509953 0.3468970 0.1857916 0.3122571
int5  0.3251571 0.3491264 0.3350408 0.1713659 0.2895193
int6  0.3049111 0.3476425 0.3224545 0.1579404 0.2691539

summary(rmdiag(corTs(ts.sim22))[btrue])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2481  0.2927  0.3381  0.3341  0.3937  0.3978 

summary(rmdiag(corTs(ts.sim1))[btrue])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3011  0.3055  0.3214  0.3218  0.3320  0.3489 

summary(rmdiag(corTs(ts.int0))[btrue])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2318  0.3543  0.3703  0.3527  0.3858  0.4212 

Global mean across interventions 0-7 and across nodes

mean(R[3:9,idx])

[1] 0.3262383

mean across interventions 0-7

colMeans(R[3:9,idx])

     1->2      2->3      3->4      1->5      4->5 
0.3769466 0.3500691 0.3521956 0.2053238 0.3466564 

Loading DGM data from Sim1 and Sim22

# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'sim1'), sprintf("Id_%03d",s), Nn)
# }
# dgm.sim1=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'sim22'), sprintf("Id_%03d",s), Nn)
# }
# dgm.sim22=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'sim22'), sprintf("Id_%03d",s), Nn, e = 26)
# }
# dgm.sim22_e26=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise01_HRF4_Mod1_Inj0_F1_60min'), sprintf("Id_%03d",s), Nn)
# }
# dgm.long=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise1_HRF4_Mod1_Inj0_F1_noise'), sprintf("Id_%03d",s), Nn)
# }
# dgm.noise=dgm.group(subj)
# 
# f=file(file.path(PATH,"results", "DGM-Sim.RData"))
# save(dgm.sim1, dgm.sim22, dgm.sim22_e26, dgm.long, dgm.noise, file = f, compress = T)
# close(f)
load(file.path(PATH, 'results', 'DGM-Sim.RData'))

Patel network analysis

set.seed(1980)
th=rand.test(ts.sim1) # get sign. thresholds
subj=list()
for (s in 1:N) {
  subj[[s]] = patel(scaleTs(ts.sim1[,,s]), TK = th$kappa, TT = th$tau)
}
patel.sim1=patel.group(subj)
th=rand.test(ts.sim22) # get sign. thresholds
subj=list()
for (s in 1:N) {
  subj[[s]] = patel(scaleTs(ts.sim22[,,s]), TK = th$kappa, TT = th$tau) # scaling is not necessary
}
patel.sim22=patel.group(subj)

spDCM

In spDCM results, rows are child nodes, columns are parent nodes

d = readMat(file.path(PATH_RES,'spDCM_Ep_A_sim1.mat'))
dcm.sim1 = d$DCM.Ep.A
d = readMat(file.path(PATH_RES,'spDCM_Ep_A_sim22.mat'))
dcm.sim22 = d$DCM.Ep.A
# s mean strenght, a thresholded adjacency
#dcm.s.sim1   = t(apply(dcm.sim1, c(1,2), mean))
#dcm.s.sim22  = t(apply(dcm.sim22, c(1,2), mean))

lingam

d = readMat(file.path(PATH_RES,'lingam_sim1.mat'))
ling.sim1 = d$LIN
d = readMat(file.path(PATH_RES,'lingam_sim22.mat'))
ling.sim22 = d$LIN
# Lingam
# as lingam only determines directionality, we supply the undirected true network
x=!btrue+t(btrue)
ling.sim1[x] = 0
ling.sim22[x] = 0

Statistical inference

stats.dgm.sim1  = binom.nettest(dgm.sim1$tam, alter = "greater", fdr = 0.05)
stats.dgm.sim22 = binom.nettest(dgm.sim22$tam, alter = "greater", fdr = 0.05)
stats.dgm.sim22_e26 = binom.nettest(dgm.sim22_e26$tam, alter = "greater", fdr = 0.05)
stats.dgm.sim22_np  = binom.nettest(dgm.sim22_e26$am, alter = "greater", fdr = 0.05)
stats.dgm.long = binom.nettest(dgm.long$tam, alter = "greater", fdr = 0.05)
stats.dgm.noise = binom.nettest(dgm.noise$tam, alter = "greater", fdr = 0.05)
# patel
stats.pat.sim1  = binom.nettest(patel.sim1$net, alter = "greater", fdr = 0.05)
stats.pat.sim22 = binom.nettest(patel.sim22$net, alter = "greater", fdr = 0.05)
# spDCM
f = 0.10
stats.DCM.sim1  = binom.nettest(dcm.sim1 > f, alter = "greater", fdr = 0.05)
stats.DCM.sim22 = binom.nettest(dcm.sim22 > f, alter = "greater", fdr = 0.05)
# Lingam
stats.ling.sim1  = binom.nettest(ling.sim1 > 0, alter = "greater", fdr = 0.05)
stats.ling.sim22 = binom.nettest(ling.sim22 > 0, alter = "greater", fdr = 0.05)

Median sensitivity and specificity

perf.dgm=list()
perf.pat=list()
perf.DCM=list()
perf.ling=list()
perf.dgm$sim1  = perf(dgm.sim1$tam, btrue)
perf.dgm$sim22 = perf(dgm.sim22$tam, btrue)
perf.dgm$long  = perf(dgm.long$tam, btrue)
perf.dgm$noise  = perf(dgm.noise$tam, btrue)
perf.dgm$sim22_e26 = perf(dgm.sim22_e26$tam, btrue)
perf.dgm$sim22_np = perf(dgm.sim22$am, btrue)
# Patel
perf.pat$sim1  = perf(patel.sim1$net, btrue)
perf.pat$sim22 = perf(patel.sim22$net, btrue)
# spDCM
perf.DCM$sim1  = perf(dcm.sim1 > f, t(btrue))
perf.DCM$sim22 = perf(dcm.sim22 > f, t(btrue))
# Lingam
perf.ling$sim1  = perf(ling.sim1 > 0, btrue)
perf.ling$sim22 = perf(ling.sim22 > 0, btrue)
table.perf=array(c(perf.dgm$sim22$tpr, perf.dgm$sim22$spc, perf.dgm$sim22$acc,
                   perf.dgm$sim1$tpr,  perf.dgm$sim1$spc,  perf.dgm$sim1$acc,
                   perf.dgm$long$tpr,  perf.dgm$long$spc,  perf.dgm$long$acc,
                   perf.dgm$noise$tpr,  perf.dgm$noise$spc,  perf.dgm$noise$acc,
                   perf.dgm$sim22_e26$tpr, perf.dgm$sim22_e26$spc, perf.dgm$sim22_e26$acc,
                   perf.dgm$sim22_np$tpr, perf.dgm$sim22_np$spc, perf.dgm$sim22_np$acc,
                   perf.pat$sim22$tpr, perf.pat$sim22$spc, perf.pat$sim22$acc,
                   perf.pat$sim1$tpr,  perf.pat$sim1$spc,  perf.pat$sim1$acc,
                   perf.DCM$sim22$tpr,  perf.DCM$sim22$spc,  perf.DCM$sim22$acc,
                   perf.DCM$sim1$tpr,  perf.DCM$sim1$spc,  perf.DCM$sim1$acc,
                   perf.ling$sim22$tpr,  perf.ling$sim22$spc,  perf.ling$sim22$acc,
                   perf.ling$sim1$tpr,  perf.ling$sim1$spc,  perf.ling$sim1$acc
                   ),
                 dim=c(3,12))
rownames(table.perf) <- c("Sensitvity", "Specificity", "Accuracy")
colnames(table.perf) <- c("DGM_Sim22", "DGM_Sim1", 'DGM_60min', 'DGM_noise',  'DGM_Sim22e26', 'DGM_Sim22np',
                          "Pat_Sim22", "Pat_Sim1", "DCM_Sim22", "DCM_Sim1", "Ling_Sim22", "Ling_Sim1")
print(table.perf, digits = 2)

            DGM_Sim22 DGM_Sim1 DGM_60min DGM_noise DGM_Sim22e26 DGM_Sim22np Pat_Sim22 Pat_Sim1 DCM_Sim22 DCM_Sim1
Sensitvity       0.70     0.50      0.91      0.70         0.67        0.84      0.42     0.47      0.36     0.58
Specificity      0.79     0.79      0.60      0.73         0.81        0.64      0.82     0.86      0.58     0.52
Accuracy         0.76     0.72      0.68      0.72         0.78        0.69      0.72     0.76      0.53     0.53
            Ling_Sim22 Ling_Sim1
Sensitvity        0.41      0.70
Specificity       0.80      0.90
Accuracy          0.70      0.85

True network detection

x=array(c(
  sum(perf.dgm$sim22$subj[,1]>=1),
  sum(perf.pat$sim22$subj[,1]>=1),
  sum(perf.ling$sim22$subj[,1]>=1),
  sum(perf.dgm$sim22$subj[,1]>=0.8),
  sum(perf.pat$sim22$subj[,1]>=0.8),
  sum(perf.ling$sim22$subj[,1]>=0.8)
  ), dim=c(3,2))/N
colnames(x)=c("5/5 nodes","4/5 nodes")
rownames(x)=c("DGM","Patel", "Lingam")
print(x)

       5/5 nodes 4/5 nodes
DGM         0.14      0.56
Patel       0.00      0.10
Lingam      0.00      0.08

Proportions

Dynamic data

# DGM
rmna(stats.dgm.sim22$adj)

     [,1] [,2] [,3] [,4] [,5]
[1,] 0.00 0.78 0.10 0.14 0.88
[2,] 0.62 0.00 0.72 0.12 0.16
[3,] 0.10 0.38 0.00 0.58 0.12
[4,] 0.08 0.12 0.38 0.00 0.52
[5,] 0.38 0.10 0.04 0.36 0.00

summary(stats.dgm.sim22$adj[btrue], na.rm = T)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.520   0.580   0.720   0.696   0.780   0.880 

# Patel
rmna(stats.pat.sim22$adj)

     [,1] [,2] [,3] [,4] [,5]
[1,] 0.00 0.52 0.42 0.12 0.66
[2,] 0.28 0.00 0.42 0.16 0.18
[3,] 0.06 0.34 0.00 0.24 0.12
[4,] 0.00 0.12 0.34 0.00 0.26
[5,] 0.16 0.18 0.06 0.16 0.00

summary(stats.pat.sim22$adj[btrue], na.rm = T)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.24    0.26    0.42    0.42    0.52    0.66 

# Lingam
rmna(stats.ling.sim22$adj)

     [,1] [,2] [,3] [,4] [,5]
[1,] 0.00 0.26  0.0  0.0 0.38
[2,] 0.74 0.00  0.4  0.0 0.00
[3,] 0.00 0.60  0.0  0.4 0.00
[4,] 0.00 0.00  0.6  0.0 0.60
[5,] 0.62 0.00  0.0  0.4 0.00

summary(stats.ling.sim22$adj[btrue], na.rm = T)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.260   0.380   0.400   0.408   0.400   0.600 

Stationary data

# DGM
rmna(stats.dgm.sim1$adj)

     [,1] [,2] [,3] [,4] [,5]
[1,] 0.00 0.38 0.06 0.16 0.56
[2,] 0.54 0.00 0.62 0.06 0.08
[3,] 0.04 0.34 0.00 0.34 0.18
[4,] 0.10 0.06 0.62 0.00 0.58
[5,] 0.38 0.04 0.06 0.40 0.00

summary(stats.dgm.sim1$adj[btrue==1], na.rm = T)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.340   0.380   0.560   0.496   0.580   0.620 

# Patel
rmna(stats.pat.sim1$adj)

     [,1] [,2] [,3] [,4] [,5]
[1,] 0.00 0.50 0.24 0.04 0.60
[2,] 0.20 0.00 0.40 0.04 0.20
[3,] 0.02 0.34 0.00 0.30 0.10
[4,] 0.04 0.10 0.38 0.00 0.54
[5,] 0.10 0.04 0.04 0.20 0.00

summary(stats.pat.sim1$adj[btrue==1], na.rm = T)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.300   0.400   0.500   0.468   0.540   0.600 

# Lingam
rmna(stats.ling.sim1$adj)

     [,1] [,2] [,3] [,4] [,5]
[1,] 0.00 0.84 0.00 0.00 0.74
[2,] 0.16 0.00 0.68 0.00 0.00
[3,] 0.00 0.32 0.00 0.68 0.00
[4,] 0.00 0.00 0.32 0.00 0.58
[5,] 0.26 0.00 0.00 0.42 0.00

summary(stats.ling.sim1$adj[btrue==1], na.rm = T)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.580   0.680   0.680   0.704   0.740   0.840 

Supplementary Table S1: time to peak

TR=0.1 # This data has a TR of 0.1
STIM_ONSET=5 # 5 sec.
table.S1=array(c(
  rowMeans(apply(ts.ssint0, c(2,3), which.max) * TR - STIM_ONSET),
  rowMeans(apply(ts.ssint1, c(2,3), which.max) * TR - STIM_ONSET),
  rowMeans(apply(ts.ssint2, c(2,3), which.max) * TR - STIM_ONSET),
  rowMeans(apply(ts.ssint3, c(2,3), which.max) * TR - STIM_ONSET),
  rowMeans(apply(ts.ssint4, c(2,3), which.max) * TR - STIM_ONSET),
  rowMeans(apply(ts.ssint5, c(2,3), which.max) * TR - STIM_ONSET),
  rowMeans(apply(ts.ssint6, c(2,3), which.max) * TR - STIM_ONSET)
  ), dim=c(Nn, 7))
table.S1sd=array(c(
  apply(apply(ts.ssint0, c(2,3), which.max) * TR - STIM_ONSET, 1, sd),
  apply(apply(ts.ssint1, c(2,3), which.max) * TR - STIM_ONSET, 1, sd),
  apply(apply(ts.ssint2, c(2,3), which.max) * TR - STIM_ONSET, 1, sd),
  apply(apply(ts.ssint3, c(2,3), which.max) * TR - STIM_ONSET, 1, sd),
  apply(apply(ts.ssint4, c(2,3), which.max) * TR - STIM_ONSET, 1, sd),
  apply(apply(ts.ssint5, c(2,3), which.max) * TR - STIM_ONSET, 1, sd),
  apply(apply(ts.ssint6, c(2,3), which.max) * TR - STIM_ONSET, 1, sd)
  ), dim=c(Nn, 7))
# Mean time to peak
print(array(as.numeric(sprintf("%.2f", table.S1)), dim=c(Nn, 7)))

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 4.30 4.57 4.79 5.06 5.29 5.51 5.71
[2,] 4.30 4.16 4.03 3.91 3.82 3.78 3.73
[3,] 4.40 4.44 4.44 4.44 4.44 4.44 4.44
[4,] 4.53 4.77 5.00 5.25 5.49 5.71 5.91
[5,] 4.32 4.09 3.98 3.88 3.80 3.76 3.73

# SD
print(array(as.numeric(sprintf("%.2f", table.S1sd)), dim=c(Nn, 7)))

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.19 0.25 0.24 0.24 0.25 0.25 0.25
[2,] 0.20 0.23 0.24 0.23 0.22 0.23 0.21
[3,] 0.28 0.23 0.23 0.23 0.23 0.23 0.23
[4,] 0.24 0.21 0.20 0.20 0.19 0.20 0.21
[5,] 0.24 0.23 0.24 0.24 0.22 0.21 0.19

Supplementary Table S2: offset relative to first simulation

print(table.S1-table.S1[,1], digits = 2)

     [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]
[1,]    0  0.276  0.494  0.758  0.996  1.210  1.416
[2,]    0 -0.138 -0.276 -0.390 -0.480 -0.520 -0.570
[3,]    0  0.046  0.046  0.046  0.046  0.046  0.046
[4,]    0  0.240  0.468  0.722  0.964  1.178  1.378
[5,]    0 -0.230 -0.346 -0.446 -0.518 -0.562 -0.588

Supplementary Table S3: total offset

M = array(NA, dim=c(3,7))
x=table.S1-table.S1[,1]
# 1 and 2
M[1,] = colSums(abs(x[c(T,T,F,F,F),]))
# 4 and 5
M[2,] = colSums(abs(x[c(F,F,F,T,T),]))
# 1 and 5
M[3,] = colSums(abs(x[c(T,F,F,F,T),]))
print(M[,2:7], digits = 3)

      [,1]  [,2] [,3] [,4] [,5] [,6]
[1,] 0.414 0.770 1.15 1.48 1.73 1.99
[2,] 0.470 0.814 1.17 1.48 1.74 1.97
[3,] 0.506 0.840 1.20 1.51 1.77 2.00

Mean cross the three edges

print(colMeans(M[,2:7]), digits = 2)

[1] 0.46 0.81 1.17 1.49 1.75 1.99

Data preparation for Figure 2

For demonstration purposes, we need to the the time series of the subject closest to the mean peak, for each node, and each intervention strength.

ix=50:110 # start is set to stimulus onset at 5s
M = array(NA, dim=c(300,Nn,N,7))
M[,,,1] = ts.ssint0
M[,,,2] = ts.ssint1
M[,,,3] = ts.ssint2
M[,,,4] = ts.ssint3
M[,,,5] = ts.ssint4
M[,,,6] = ts.ssint5
M[,,,7] = ts.ssint6
S=array(NA, dim=c(Nn,7))
for (n in 1:Nn){
  for (i in 1:7){
  S[n,i]=which.min(abs(table.S1[n,i]-apply(M[ix,n,,i], 2, which.max)*0.1))
  }
}
n=5
for (i in 1:7) {
  dt = abs(table.S1[n,i]-apply(M[ix,n,,i], 2, which.max)*0.1)
  m=min(dt)
  #print(which(m==dt))
}
# replace some subjects with others that have same time to peak 
S[1,3]=10
S[2,c(4,6)]=c(12,19)
S[3,2:7]=c(15,16,20,28,29,31)
S[4,7]=6
S[5,c(1,2,7)]=c(9,31,6)

Figure 2: Interventions

img = readPNG(file.path(PATH_FIG, "fig-interventions-page001.png"))
g = rasterGrob(img, interpolate=T)
pA = ggplot() + annotation_custom(g) + theme(plot.title = element_text(size=12)) +
  ggtitle('HRF lag intervention')
pB=gplotMat(R[,idx], lim = c(0.1, 0.5), colMapLabel = expression("Pearson\'s"~italic(r)), barWidth = 0.2,
            title = "node correlations", titleTextSize = 12) + xlab("Node pairs") +
  ylab("Dataset") + scale_x_discrete(limits=c("1\n2","2\n3","3\n4", "4\n5", "1\n5")) +
  scale_y_discrete(limits=c("Sim1", "Sim22", "60 min.","< 0.4 s", "0.4 s", "0.8 s",
                            "1.1 s", "1.4 s", "1.7 s", "1.9 s")) +
  theme(axis.text.y = element_text(size=11))
l=length(ix)
offset=as.factor(c(rep("< 0.4 s",l), rep("0.4 s",l), rep("0.8 s",l), rep("1.1 s",l),
                   rep("1.4 s",l), rep("1.7 s",l), rep("1.9 s",l)))
p=list()
mylegend = c(rep("none",4), "right")
mytitles = c("node 1 +delay", "node 2 -delay", "node 3", "node 4 +delay", "node 5 -delay")
m = array(NA, dim=c(l,7))
for (n in 1:Nn){ 
  x=array(NA, dim=c(l,7))
  for (i in 1:7) {
    x[,i] = M[ix,n,S[n,i],i]
    m[,i] = rep(table.S1[n,i],l)
  }
  x=melt(x)
  x$offset=offset
  x$m=c(m)
  
  p[[n]] = ggplot(x, aes(x=Var1*TR, y=value, group=Var2, colour=offset)) + geom_line(size=1) + 
    ggtitle(mytitles[n]) + xlab("time (s)") +
    theme(legend.position=mylegend[n], plot.title = element_text(size=12)) + 
    geom_vline(data = x, aes(xintercept = m, color=offset))
}
top=plot_grid(pA, pB, labels=c("A", "B"), ncol = 2, nrow = 1, rel_widths = c(0.8, 1))
mid=plot_grid(p[[1]], p[[2]], p[[3]], labels="C", ncol = 3, nrow = 1, rel_widths = c(1, 1, 1))
bot=plot_grid(p[[4]], p[[5]], ncol = 2, nrow = 1, rel_widths = c(0.7, 1))
plot_grid(top, mid, bot, ncol=1, nrow=3, rel_heights = c(1, 0.7, 0.7))

ggsave(path = PATH_FIG, "Fig2.png")

Figure: Various delays at node 5

# ts x nodes x subj
s=1
n=4
idx=50:180
x = cbind(ts.ssint0[idx,n,s], ts.ssint1[idx,n,s], ts.ssint2[idx,n,s],
          ts.ssint3[idx,n,s], ts.ssint4[idx,n,s], ts.ssint5[idx,n,s],
          ts.ssint6[idx,n,s])
#plot.ts(x)
ggplot(melt(x), aes(x=Var1/10, y=value, group=Var2, colour=Var2)) +
  geom_line(size=1) + ggtitle(mytitles[n]) + xlab("time (s)") 

Figure: Estimates of theta as function of time

# example dataset and node 2 has node 1 as parent 1->2
s=10
node=2
pars=1
Nt=nrow(ts.sim22)
TR=3
# dgm.sim22$tam[,,s]
Ft=array(1,dim=c(Nt,length(pars)+1))
Ft[,2:ncol(Ft)]=ts.sim22[,pars,s] # selects parents
Yt=ts.sim22[,node,s]
# get df corresponding to parent model pars
df = getModel(dgm.sim22$models[,,node,s], pars)[Nn+2] 
fit=dlm.lpl(Yt, t(Ft), delta = df)
y = dlm.retro(fit$mt, fit$CSt, fit$RSt, fit$nt, fit$dt)
bold=ts.sim22[,c(1,2),s]
theta=cbind(y$smt[2,])
 p1 = ggplot(melt(bold), aes(x = Var1*TR, y = value, group=Var2, color=as.factor(Var2))) + geom_line() +
   theme_minimal() + ggtitle("simulated fMRI time series of two nodes") + 
   scale_color_discrete(name = "node") + xlab("seconds") + ylab("")
 p2 = ggplot(melt(theta), aes(x = Var1*TR, y = value, group=Var2, color=as.factor("1"))) + geom_line() +
   theme_minimal() + ggtitle("connectivity thrength over time") + scale_color_discrete(name = "theta") +
   xlab("seconds") + ylab("theta")
 
plot_grid(p1, p2, ncol = 1, nrow = 2, rel_widths = c(1, 1))

ggsave(path = PATH_FIG, "Theta.png")

Saving 5 x 3 in image

Figure 5: DGM vs. other methods

s = 0.2 # spacing
pA=ggplot(melt(table.perf[1:3,c(1,5,6,7,11)]), aes(x=Var2, y=value, fill=Var1)) +
  scale_x_discrete(labels=c("DGM_Sim22" = "DGM\ne=20", "DGM_Sim22e26" = "DGM\ne=26",
                            "DGM_Sim22np" = "DGM\ne=0", "Pat_Sim22" = "Patel",
                            "DCM_Sim22" = "spDCM", "Ling_Sim22" = "Ling")) +
  geom_bar(stat="identity", position=position_dodge()) + guides(fill=FALSE) + ylab("Proportion") +
  theme(axis.title.x=element_blank(), axis.text.x = element_text(size=10), axis.text.y = element_text(size=9),
        axis.title.y = element_text(size=10)) + 
  ggtitle("dynamic nodes") + coord_cartesian(ylim=c(0.4,1)) + scale_fill_brewer(palette="Set1")
pB=ggplot(melt(table.perf[1:3,c(2,8,12)]), aes(x=Var2, y=value, fill=Var1)) +
  scale_x_discrete(labels=c("DGM_Sim1" = "DGM\ne=20", "Pat_Sim1" = "Patel",
                            "DCM_Sim1" = "spDCM", "Ling_Sim1" = "Ling")) + 
  geom_bar(stat="identity", position=position_dodge()) + ylab("Proportion") + ggtitle("stationary nodes") +
  coord_cartesian(ylim=c(0.4,1)) + guides(fill=guide_legend(title="")) +
  theme(axis.title.x=element_blank(), axis.text.x = element_text(size=10), axis.text.y = element_text(size=9),
        axis.title.y = element_text(size=10), legend.text=element_text(size=10)) +
  scale_fill_brewer(palette="Set1")
pTop = plot_grid(pA, pB,  nrow=1, ncol=2, rel_widths = c(0.85, 1), labels = c("A", "B")) +
  theme(plot.margin = unit(c(s, 0, s, 0), "cm"))
pC1 = gplotMat(stats.dgm.sim22$adj, 'DGM', '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pC2 = gplotMat(stats.pat.sim22$adj, 'Patel', '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
#pC3 = gplotMat(rmdiag(t(stats.DCM.sim22$adj)), "spDCM", '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pC4 = gplotMat(rmdiag(stats.ling.sim22$adj), "Lingam", '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pC5 = gplotMat(stats.dgm.sim22$adj_fdr, 'DGM (FDR)', '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pC6 = gplotMat(stats.pat.sim22$adj_fdr, 'Patel (FDR)', '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
#pC7 = gplotMat(rmdiag(t(stats.DCM.sim22$adj_fdr)), "spDCM (FDR)", '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pC8 = gplotMat(rmdiag(stats.ling.sim22$adj_fdr), "Lingam (FDR)", '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pD1 = gplotMat(stats.dgm.sim1$adj, 'DGM', '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pD2 = gplotMat(stats.pat.sim1$adj, 'Patel', '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
#pD3 = gplotMat(rmdiag(t(stats.DCM.sim1$adj)), "spDCM", '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pD4 = gplotMat(rmdiag(stats.ling.sim1$adj), "Lingam", '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pD5 = gplotMat(stats.dgm.sim1$adj_fdr, 'DGM (FDR)', '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pD6 = gplotMat(stats.pat.sim1$adj_fdr, 'Patel (FDR)', '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pD7 = gplotMat(rmdiag(t(stats.DCM.sim1$adj_fdr)), "spDCM (FDR)", '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pD8 = gplotMat(rmdiag(stats.ling.sim1$adj_fdr), "Lingam (FDR)", '%', titleTextSize = 10, axisTextSize=10, textSize = 10, barWidth = 0.2)
pMid  = plot_grid(pC1, pC2, pC4, pC5, pC6, pC8, nrow=2, ncol=3, rel_widths = c(1,1,1)) +
  theme(plot.margin = unit(c(s, 0, s, 0), "cm"))
pBott = plot_grid(pD1, pD2, pD4, pD5, pD6, pD8, nrow=2, ncol=3, rel_widths = c(1,1,1)) +
  theme(plot.margin = unit(c(s, 0, s, 0), "cm"))
plot_grid(pTop, pMid, pBott, ncol = 1, nrow = 3, rel_heights = c(0.55,1,1),
          labels = c("", "C dynamic nodes", "D stationary nodes"),
          vjust = 0.6, hjust = -0.1)

ggsave(path = PATH_FIG, "Fig5.png")

Saving 5.5 x 7.79 in image

Common network

Maximizing LPLs across datasets

# dim(dgm.sim22$models)
# 1. sum all LPLs across subjects
# 2. for each child node, maximize across models
idx = apply(apply(dgm.sim22$models[Nn+1,,,], c(1,2), sum), 2, which.max)
# create network matrix
M = array(0, dim=c(Nn, Nn))
for (n in 1:Nn) {
  M[dgm.sim22$models[2:Nn,idx[n],n,1], n] = 1
}
gplotMat(M, title = "Common net", hasColMap = F)

Loading DGM of 7 HRF datasets

# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise01_HRF1_Mod1_Inj0_F1'),
#                            sprintf("Id_%03d",s), Nn)
#   subj[[s]]$thr = pruning(subj[[s]]$adj, subj[[s]]$models, winner = subj[[s]]$winner, e = 20)
# }
# dgm.int0=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F13'),
#                            sprintf("Id_%03d",s), Nn)
# }
# dgm.int1=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F16'),
#                            sprintf("Id_%03d",s), Nn)
# }
# dgm.int2=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F20'),
#                            sprintf("Id_%03d",s), Nn)
# }
# dgm.int3=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F24'),
#                            sprintf("Id_%03d",s), Nn)
# }
# dgm.int4=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F28'),
#                            sprintf("Id_%03d",s), Nn)
# }
# dgm.int5=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_NET,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F32'),
#                            sprintf("Id_%03d",s), Nn)
# }
# dgm.int6=dgm.group(subj)
# 
# f=file(file.path(PATH,"results", "DGM-Sim_hrf.RData"))
# save(dgm.int0, dgm.int1, dgm.int2, dgm.int3, dgm.int4, dgm.int5, dgm.int6, file = f, compress = T)
# close(f)
load(file.path(PATH, 'results', 'DGM-Sim_hrf.RData'))

Investigate discount factor

node=as.factor(c(rep(1,N),rep(2,N),rep(3,N),rep(4,N),rep(5,N)))
d = list()
d[[1]]=data.frame(df=c(dgm.sim1$df_),  node=node)
d[[2]]=data.frame(df=c(dgm.sim22$df_), node=node)
d[[3]]=data.frame(df=c(dgm.int0$df_),  node=node)
d[[4]]=data.frame(df=c(dgm.int1$df_),  node=node)
d[[5]]=data.frame(df=c(dgm.int2$df_),  node=node)
d[[6]]=data.frame(df=c(dgm.int3$df_),  node=node)
d[[7]]=data.frame(df=c(dgm.int4$df_),  node=node)
d[[8]]=data.frame(df=c(dgm.int5$df_),  node=node)
d[[9]]=data.frame(df=c(dgm.int6$df_),  node=node)
p = list()
for (i in 1:9) {
  p[[i]] = ggplot(d[[i]], aes(x=node, y=df)) + geom_boxplot(width=0.4) + ggtitle(str_int[i]) +
    geom_point(shape=1, color="gray70", size=0.5, position = position_jitter(width = 0.2, height = 0.0))
}
  
plot_grid(plotlist = p, ncol = 3, nrow = 3)

DFs with parents only

x = t(apply(dgm.sim22$tam, 3, colSums)) # Subj x no. of parents
df.22 = dgm.sim22$df_
df.22[x == 0] = NA
summary(colMeans(df.22, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.7991  0.8411  0.8641  0.8598  0.8934  0.9012 

x = t(apply(dgm.sim1$tam, 3, colSums)) # Subj x no. of parents
df.1 = dgm.sim1$df_
df.1[x == 0] = NA
summary(colMeans(df.1, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.9928  0.9930  0.9940  0.9947  0.9955  0.9981 

x = t(apply(dgm.int0$am, 3, colSums)) # Subj x no. of parents
df.0 = dgm.int0$df_
df.0[x == 0] = NA
summary(colMeans(df.0, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.6492  0.6572  0.6664  0.6688  0.6731  0.6981 

x = t(apply(dgm.int6$am, 3, colSums)) # Subj x no. of parents
df.6 = dgm.int6$df_
df.6[x == 0] = NA
summary(colMeans(df.6, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.5769  0.6453  0.7080  0.6894  0.7518  0.7651 

  p1 = ggplot(melt(df.22), aes(x=Var2, y=value, group=Var2)) + geom_boxplot(width=0.4) +
    geom_point(shape=1, color="gray70", size=0.5, position = position_jitter(width = 0.2, height = 0.0))
  p2 = ggplot(melt(df.0), aes(x=Var2, y=value, group=Var2)) + geom_boxplot(width=0.4) +
    geom_point(shape=1, color="gray70", size=0.5, position = position_jitter(width = 0.2, height = 0.0))
  
plot_grid(p1, p2, ncol = 2, nrow = 1)

Stats for DGM 7 HRF datasets

stats.dgm.int0 = binom.nettest(dgm.int0$tam, alter = "greater", fdr = 0.05)
stats.dgm.int1 = binom.nettest(dgm.int1$tam, alter = "greater", fdr = 0.05)
stats.dgm.int2 = binom.nettest(dgm.int2$tam, alter = "greater", fdr = 0.05)
stats.dgm.int3 = binom.nettest(dgm.int3$tam, alter = "greater", fdr = 0.05)
stats.dgm.int4 = binom.nettest(dgm.int4$tam, alter = "greater", fdr = 0.05)
stats.dgm.int5 = binom.nettest(dgm.int5$tam, alter = "greater", fdr = 0.05)
stats.dgm.int6 = binom.nettest(dgm.int6$tam, alter = "greater", fdr = 0.05)

Sensitivity and specificity of DGM in 7 HRF datasets

perf.dgm$int0 = perf(dgm.int0$tam, btrue)
perf.dgm$int1 = perf(dgm.int1$tam, btrue)
perf.dgm$int2 = perf(dgm.int2$tam, btrue)
perf.dgm$int3 = perf(dgm.int3$tam, btrue)
perf.dgm$int4 = perf(dgm.int4$tam, btrue)
perf.dgm$int5 = perf(dgm.int5$tam, btrue)
perf.dgm$int6 = perf(dgm.int6$tam, btrue)

Estimate null sensitivity and specificity

R=array(NA, dim=c(Nn, Nn, N))
# correlation matrix for each data set
for (s in 1:N) {
  R[,,s]=cor(ts.int0[,,s])
}
summary(R[btrue])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.6981  0.1872  0.3735  0.3527  0.5689  0.8418 

summary(R[rmna(btrue) + t(rmna(btrue))==0])

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.46008  0.06037  0.29945  0.43341  1.00000  1.00000 

am = R > 0.30
# threshold these correlation matrices from null network, or
# true edges are know, but not the direction
# am = array(rep(rmna(btrue) + t(rmna(btrue)), N), dim=c(Nn,Nn,N))
#xfalse_ = array(rep(rmna(btrue) + t(rmna(btrue)), N), dim=c(Nn,Nn,N)) == 0
#summary(R[xfalse_])
#R_=R > 0.2
for (s in 1:N) {
  for (i in 1:Nn) {
    for (j in 1:Nn) {
      if (i != j & i > j) {
        x=sample(1:3, 1, replace = T) # case 3 is retaining both edges (unidrected model)
        if (x == 1) {
          am[i,j,s] = 0 # remove false edge
        } else if (x == 2) {
          am[j,i,s] = 0 # remove true edge
        }
      }
    }
  }
}
perf.dgm$null = perf(am, btrue)

c-sensitivity and d-accuracy

c-sensitivity DGM vs Patel

This compares overall c-sensitivity of DGM vs Patel with stationary (sim1) and time-varying (sim22) data. TR is 3 s.

# for c-sensitivity direction is irrelevant so we use the symmetric function.
x = array(c(sum(symmetric(dgm.sim1$tam)[btrue])/(Nn*N),
            sum(symmetric(dgm.sim22$tam)[btrue])/(Nn*N),
            sum(patel.sim1$tkappa[btrue]>0)/(Nn*N),
            sum(patel.sim22$tkappa[btrue]>0)/(Nn*N)),
          dim = c(2,2))
colnames(x) =  c("DGM", "Patel")
rownames(x) =  c("Sim1", "Sim22")
print(x)

        DGM Patel
Sim1  0.912 0.740
Sim22 0.888 0.716

c-sensitivity for 60 min. run

With time-varying data

sum(symmetric(dgm.long$tam)[btrue])/(Nn*N)

[1] 0.98

c-sensitivity for DGM varying HRF responses

This compares overall c-sensitivity of 7 different interventions to lag the HRF response.

# for c-sensitivity direction is irrelevant so we extract both the true network and the transposed network of the opposite direction and do a logical or (max).
x = array(c(sum(symmetric(dgm.int0$tam)[btrue])/(Nn*N),
            sum(symmetric(dgm.int1$tam)[btrue])/(Nn*N),
            sum(symmetric(dgm.int2$tam)[btrue])/(Nn*N),
            sum(symmetric(dgm.int3$tam)[btrue])/(Nn*N),
            sum(symmetric(dgm.int4$tam)[btrue])/(Nn*N),
            sum(symmetric(dgm.int5$tam)[btrue])/(Nn*N),
            sum(symmetric(dgm.int6$tam)[btrue])/(Nn*N)),
          dim = c(1,7))
colnames(x) =  c("int0", "int1", "int2", "int3", "int4", "int5", "int6")
print(x)

     int0  int1  int2  int3  int4  int5  int6
[1,] 0.92 0.944 0.912 0.876 0.836 0.768 0.708

d-accuracy for DGM vs Patel

x = array(c(sum(dgm.sim1$tam[btrue])/(Nn*N),
            sum(dgm.sim22$tam[btrue])/(Nn*N),
            sum(patel.sim1$net[btrue])/(Nn*N),
            sum(patel.sim22$net[btrue])/(Nn*N),
            sum(patel.sim1$tau[btrue]>0)/(Nn*N),
            sum(patel.sim22$tau[btrue]>0)/(Nn*N),
            sum(am[btrue])/(Nn*N),
            sum(am[btrue])/(Nn*N)),
          dim=c(2,4))
colnames(x) <- c("DGM", "Patel (sign. kappa and tau)", "Patel tau", "null")
rownames(x) <- c("Sim1", "Sim22")
print(x, digits = 3)

        DGM Patel (sign. kappa and tau) Patel tau  null
Sim1  0.496                       0.468     0.344 0.384
Sim22 0.696                       0.420     0.332 0.384

d-accuracy for long 60min. simulation

sum(dgm.long$tam[btrue])/(Nn*N)

[1] 0.912

d-accuracy for DGM varying HRF responses

x = array(c(sum(dgm.int0$tam[btrue])/(Nn*N),
            sum(dgm.int1$tam[btrue])/(Nn*N),
            sum(dgm.int2$tam[btrue])/(Nn*N),
            sum(dgm.int3$tam[btrue])/(Nn*N),
            sum(dgm.int4$tam[btrue])/(Nn*N),
            sum(dgm.int5$tam[btrue])/(Nn*N),
            sum(dgm.int6$tam[btrue])/(Nn*N)),
          dim = c(1,7))
colnames(x) =  c("int0", "int1", "int2", "int3", "int4", "int5", "int6")
print(x)

      int0  int1  int2  int3  int4  int5  int6
[1,] 0.796 0.772 0.724 0.684 0.616 0.552 0.484

Median Sensitivity and Specificity for HRF Interventions

res = array(NA, dim=c(3,7))
res[1,] = c(perf.dgm$int0$tpr, perf.dgm$int1$tpr, perf.dgm$int2$tpr, perf.dgm$int3$tpr,
            perf.dgm$int4$tpr, perf.dgm$int5$tpr, perf.dgm$int6$tpr)
res[2,] = c(perf.dgm$int0$spc, perf.dgm$int1$spc, perf.dgm$int2$spc, perf.dgm$int3$spc,
            perf.dgm$int4$tpr, perf.dgm$int5$spc, perf.dgm$int6$spc)
res[3,] = c(perf.dgm$int0$acc, perf.dgm$int1$acc, perf.dgm$int2$acc, perf.dgm$int3$acc,
            perf.dgm$int4$acc, perf.dgm$int5$acc, perf.dgm$int6$acc)
colnames(res) <- c("<0.4s", "0.4s", "0.8s", "1.1s", "1.4s", "1.7s", "1.9s")
rownames(res) <- c("Sensitivity", "Specificity", "Accuracy")
print(res)

                <0.4s      0.4s      0.8s      1.1s  1.4s  1.7s  1.9s
Sensitivity 0.7960000 0.7720000 0.7240000 0.6840000 0.616 0.552 0.484
Specificity 0.6893333 0.6666667 0.6666667 0.6653333 0.616 0.652 0.652
Accuracy    0.7160000 0.6930000 0.6810000 0.6700000 0.646 0.627 0.610

summary(t(res))

  Sensitivity      Specificity        Accuracy     
 Min.   :0.4840   Min.   :0.6160   Min.   :0.6100  
 1st Qu.:0.5840   1st Qu.:0.6520   1st Qu.:0.6365  
 Median :0.6840   Median :0.6653   Median :0.6700  
 Mean   :0.6611   Mean   :0.6583   Mean   :0.6633  
 3rd Qu.:0.7480   3rd Qu.:0.6667   3rd Qu.:0.6870  
 Max.   :0.7960   Max.   :0.6893   Max.   :0.7160  

Figure 6: DGM sensitivity and specificity for the 7 HRF datasets

ggplot(melt(res), aes(x=Var2, y=value, group=Var1, color=Var1)) + 
  geom_point(size=2) + geom_line(size=0.5) + ylim(c(0,0.8)) +
  theme(axis.text.x = element_text(size=9), panel.grid.major = element_line(colour = "gray70", linetype = "dotted")) + 
  guides(color=guide_legend(title="")) +
  xlab("Intervention")

ggsave(path = PATH_FIG, "Fig6.png")

Saving 4.2 x 2 in image

With a 3T scanner thermal noise is below 1%, usually ~0.2%

p1=gplotMat(stats.dgm.sim1$adj,  title = "sim1")
p2=gplotMat(stats.dgm.sim22$adj, title = "sim22")
p3=gplotMat(stats.dgm.long$adj,  title = "60 min.")
p4=gplotMat(stats.dgm.noise$adj, title = "0.3% noise")
plot_grid(p1, p2, p3, p4, ncol=2, nrow=2)

Sensitivity and Specificity for example networks

a = array(0, dim=c(Nn,Nn))
a[1,2] = a[2,3] = a[3,4] = a[1,5] = a[2,1] = 1
pa = perf(a, btrue)
p1=gplotMat(a, title = "a")
b = array(0, dim=c(Nn,Nn))
b[1,2] = b[2,3] = b[3,4] = b[4,5] = b[1,5] = 1
b[2,1] = b[3,2] = b[4,3] = b[5,4] = b[5,1] = 1
pb = perf(b, btrue)
p2=gplotMat(b, title = "b")
c = array(0, dim=c(Nn,Nn))
c[1,2] = c[2,1] = c[2,3] = c[3,2] = c[5,1] = 1
pc = perf(c, btrue)
p3=gplotMat(c, title = "c")
p = perf(btrue, btrue)
p0=gplotMat(atrue, title = "true")
plot_grid(p0, p1, p2, p3, ncol=2, nrow=2)

result = rbind(p$subj, pa$subj, pb$subj, pc$subj)
rownames(result) = c("true", "a", "b", "c")
print(round(result, digits = 2))

     tpr  spc ppv  npv  fpr fnr fdr  acc
true 1.0 1.00 1.0 1.00 0.00 0.0 0.0 1.00
a    0.8 0.93 0.8 0.93 0.07 0.2 0.2 0.90
b    1.0 0.67 0.5 1.00 0.33 0.0 0.5 0.75
c    0.4 0.80 0.4 0.80 0.20 0.6 0.6 0.70

Neural lag 50 ms and 500 ms

res=200
stim = 8401/res  # onset simimulus
z = 8410/res # 66% of amplitude z output with 50 ms
z2 = 8507/res # 66% of amplitude z output with 500 ms
print(z-stim)

[1] 0.045

print(z2-stim)

[1] 0.53

Figure: pruning example

n=10
s = read.subject(file.path(PATH_NET,'sim22'), sprintf("Id_%03d",n), Nn)
o0  = pruning(s$adj, s$models, winner = s$winner, e = 0)
o5  = pruning(s$adj, s$models, winner = s$winner, e = 5)
o10 = pruning(s$adj, s$models, winner = s$winner, e = 10)
o20 = pruning(s$adj, s$models, winner = s$winner, e = 20)
p1= gplotMat(o0$am, hasColMap = F, title = "e=0", titleTextSize = 10)
p2= gplotMat(o5$am, hasColMap = F, title = "e=5", titleTextSize = 10)
p3= gplotMat(o10$am, hasColMap = F, title = "e=10", titleTextSize = 10)
p4= gplotMat(o20$am, hasColMap = F, title = "e=20", titleTextSize = 10)
plot_grid(p1, p2, p3, p4, ncol=4, nrow = 1)

ggsave(path = PATH_FIG, "PruningExample.png")

Saving 6 x 1.7 in image