DGM-Stepwise

Install required packages

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

Load libraries

library(multdyn)
library(R.matlab)
library(testit)
library(ggplot2)
library(cowplot)
library(reshape2)
library(png)

Variables

PATH_HOME = "/home/simon"
PATH = file.path(PATH_HOME, "Dropbox", "Data", "DGM")  # Project path
PATH_FIG  = file.path(PATH, 'figures') # path where figures 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

5-node network simulation

Estimate networks

N=50 # Number of simulated subjects/datasets
Nn=5 # Number of nodes
# ts.int1 = readMat(file.path(PATH_TS,'Nn5_TR2_Noise01_HRF1_Mod1_Inj1_F13.mat'))$gfy2s
# ex = fw = bw = bo = list()
# for (i in 1:N) {
#   print(sprintf("Estimating subject %03d", i))
#   ex[[i]] = subject(ts.int1[,,i], method = "exhaustive")
#   fw[[i]] = subject(ts.int1[,,i], method = "forward")
#   bw[[i]] = subject(ts.int1[,,i], method = "backward")
#   bo[[i]] = subject(ts.int1[,,i], method = "both")
# }
# 
# f=file(file.path(PATH,"results", "DGM-Step.RData"))
# save(ex, fw, bw, bo, file = f, compress = T)
# close(f)
load(file.path(PATH, 'results', 'DGM-Step.RData'))

Compare stepwise

g.ex = dgm.group(ex)
g.fw = dgm.group(fw)
g.bw = dgm.group(bw)
g.bo = dgm.group(bo)
prc = array(NA,dim = c(1,3))
colnames(prc) = c("forward", "backward", "both")
prc[1] = (sum(g.ex$am==g.fw$am)-(Nn*N))/(Nn*(Nn-1)*N)
prc[2] = (sum(g.ex$am==g.bw$am)-(Nn*N))/(Nn*(Nn-1)*N)
prc[3] = (sum(g.ex$am==g.bo$am)-(Nn*N))/(Nn*(Nn-1)*N)
prc=round(prc*100, 2)
print(prc)

     forward backward both
[1,]    97.9     97.8 99.8

Stepwise both forward and backward combined is best and almost identical to exhaustive search.

10-node network ICA fMRI

Estimating stepwise 10 nodes

# load(file.path(PATH, 'results', 'RSN10-ts.RData'))
# subj_r1 = vector(mode = "list", length = N)
# subj_r2 = vector(mode = "list", length = N)
# subj_r3 = vector(mode = "list", length = N)
# subj_r4 = vector(mode = "list", length = N)
# 
# for (s in 1:N) {
#   subj_r1[[s]] = read.subject(PATH_NET_BO, sprintf("%s_Run_%03d", SUBJECTS[s], 1), Nn)
#   subj_r2[[s]] = read.subject(PATH_NET_BO, sprintf("%s_Run_%03d", SUBJECTS[s], 2), Nn)
#   subj_r3[[s]] = read.subject(PATH_NET_BO, sprintf("%s_Run_%03d", SUBJECTS[s], 3), Nn)
#   subj_r4[[s]] = read.subject(PATH_NET_BO, sprintf("%s_Run_%03d", SUBJECTS[s], 4), Nn)
#   print(sprintf("Subject %03d loaded",s))
# }
# 
# dgm.net10.bo.r1 = dgm.group(subj_r1)
# dgm.net10.bo.r2 = dgm.group(subj_r2)
# dgm.net10.bo.r3 = dgm.group(subj_r3)
# dgm.net10.bo.r4 = dgm.group(subj_r4)
# 
# f=file(file.path(PATH,"results", "DGM-RSN10_bo.RData"))
# save(dgm.net10.bo.r1, dgm.net10.bo.r2, dgm.net10.bo.r3, dgm.net10.bo.r4, file = f, compress = T)
# close(f)
load(file.path(PATH, 'results', 'DGM-RSN10_bo.RData'))
load(file.path(PATH, 'results', 'DGM-RSN10.RData'))

Compare stepwise

N=500
Nn=10
prc10 = array(NA,dim = c(1,4))
colnames(prc10) = c("run 1", "run 2", "run 3", "run 4")
prc10[1] = (sum(dgm.net10.r1$am==dgm.net10.bo.r1$am)-(Nn*N))/(Nn*(Nn-1)*N)
prc10[2] = (sum(dgm.net10.r2$am==dgm.net10.bo.r2$am)-(Nn*N))/(Nn*(Nn-1)*N)
prc10[3] = (sum(dgm.net10.r3$am==dgm.net10.bo.r3$am)-(Nn*N))/(Nn*(Nn-1)*N)
prc10[4] = (sum(dgm.net10.r4$am==dgm.net10.bo.r4$am)-(Nn*N))/(Nn*(Nn-1)*N)
prc10=round(prc10*100, 2)
print(prc10)

     run 1 run 2 run 3 run 4
[1,] 99.36 99.52 99.49 99.52

Plot

p1=ggplot(melt(prc), aes(x=Var2, y=value, fill=Var2)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=as.vector(prc)), size=4, vjust=-0.3) +
  ggtitle("5-node network") + ylim(c(0,105)) + ylab("percent") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        legend.title=element_blank())
p2=ggplot(melt(prc10), aes(x=Var2, y=value, fill=Var2)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=as.vector(prc10)), size=4, vjust=-0.3) +
  ggtitle("10-node RSNs") + ylim(c(0,105)) + ylab("percent") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        legend.title=element_blank())
plot_grid(p1, p2, ncol=2, rel_widths = c(1, 1))