Model

For details on the model, please see SIEM model.

Real Data Experiments

In this notebook, we determine whether there exists a difference in connectivity \(p_{homo}\) homotopically (same region, opposite hemisphere) vs \(p_{hetero}\) heterotopically (different region) connectivity within a particular modality. We consider \(\delta_{x} = p_{homo} - p_{hetero}\) to be the difference in connectivity for a graph from a graph or collection of graphs from a particular modality \(x\).

Test

Our test for this notebook is as follows:

\[\begin{align*} H_0: \delta_d \leq \delta_f \\ H_A: \delta_d > \delta_f \end{align*}\]

in words, whether the connectivity homotopically exceeds the connectivity heterotopically. We will do this in 2 ways:

  • An \(n\) sample test to determine whether the difference in connectivity homotopically vs heterotopically in the functional connectomes exceeds that of the diffusion connectomes at the population level.
  • A \(2\)-sample test to determine whether the difference in connectivity homotopically vs heterotopically in the functional connectomes exceeds that of thet diffusion connectomes at the individual graph level (1 diffusion and 1 functional graph from the same subject).

We will perform this experiment for both dMRI and fMRI-derived connectomes.

Raw Data

For the data, we compute the weighted mean functional (rank of each edge) and diffusion (number of fibers). For the functional connectome, we threshold such that the largest 50% of edges are set to connected, and the smallest 50% set to disconnected. For the diffusion (which are natively sparse) we just threshold edges that are present to connected, and edges that are not present to disconnected (threshold about 0).

The data below can be downloaded and moved to appropriate folders as follows (note that the below section requires sudo access):

sudo mkdir /data/
sudo chmod -R 777 /data

cd /data
wget http://openconnecto.me/mrdata/share/derivatives/dwi_edgelists.tar.gz
wget http://openconnecto.me/mrdata/share/derivatives/fmri_edgelists.tar.gz
wget http://openconnecto.me/mrdata/share/connectome_stats/connectome_stats.zip

mkdir -p /data/connectome_stats /data/all_mr /data/all_mr/dwi/edgelists /data/all_mr/fmri/ranked/edgelists
mv dwi_edgelists.tar.gz /data/dwi/edgelists
cd /data/dwi/edgelists
tar -xvzf dwi_edgelists.tar.gz
mv /data/fmri_edgelists.tar.gz /data/fmri/ranked/edgelists
cd /data/fmri/ranked/edgelists
tar -xvzf fmri_edgelists.tar.gz
mv /data/connectome_stats.zip /data/connectome_stats.zip
cd /data/connectome_stats
unzip connectome_stats.zip
basepath = '/data/connectome_stats/'
fmri_gr = read_graph(file.path(basepath, 'fmrimean_1709.edgelist'), format="ncol")
vset <- V(fmri_gr)
ordered_v <- order(vset)
fmri_gr = read_graph(file.path(basepath, 'fmrimean_1709.edgelist'), format="ncol", predef=ordered_v)
fmri_mean = get.adjacency(fmri_gr, type="both", sparse=FALSE, attr='weight')
dwi_gr = read_graph(file.path(basepath, 'dwimean_2861.edgelist'), format="ncol", predef=ordered_v)
dwi_mean = get.adjacency(dwi_gr, type="both", sparse=FALSE, attr='weight')

fmri_thresh = thresh_matrix(fmri_mean)
dwi_thresh = thresh_matrix(dwi_mean, thresh=0)

gs.plot.plot_matrix(fmri_thresh, title = "Mean Thresholded Functional Connectome", legend.name = "connection", ffactor = TRUE) +
  theme_bw()

gs.plot.plot_matrix(dwi_thresh, title = "Mean Thresholded DWI Connectome", legend.name = "connection", ffactor = TRUE) +
  theme_bw()

Blocked Data

here, we will compute the probability of an edge existing in each of 4 quadrants (2 homotopic quadrants; 2 heterotopic quadrants):

group1 = 1:35
group2 = 36:70
groups = list(group1, group2)
fmri_block = block_data(fmri_thresh, groups)
dwi_block = block_data(dwi_thresh, groups)

fmriu.plot.plot_graph(fmri_block, title = "Blocked Functional Connectome", xlabel = "Hemisphere",
                      ylabel="Hemisphere", include_diag = TRUE, legend.name = "p")

fmriu.plot.plot_graph(dwi_block, title = "Blocked DWI Connectome", xlabel = "Hemisphere",
                      ylabel="Hemisphere", include_diag = TRUE, legend.name = "p")

Diffusion

nroi <- 70
dwi.dsets = c('BNU1', 'BNU3', 'HNU1', 'KKI2009', 'NKI1', 'NKIENH', 'MRN1313', 'Templeton114', 'Templeton255', 'SWU4')
dwi.atlas = 'desikan'
dwi.basepath = '/data/all_mr/dwi/edgelists'

graphobj = fmriu.io.collection.open_graphs(basepath = dwi.basepath, atlases = dwi.atlas, datasets = dwi.dsets,
                                           gname = 'graphs', fmt='edgelist', rtype = 'array')
## [1] "opening graphs for BNU1 dataset and desikan parcellation atlas..."
## [1] "opening graphs for BNU3 dataset and desikan parcellation atlas..."
## [1] "opening graphs for HNU1 dataset and desikan parcellation atlas..."
## [1] "opening graphs for KKI2009 dataset and desikan parcellation atlas..."
## [1] "opening graphs for NKI1 dataset and desikan parcellation atlas..."
## [1] "opening graphs for NKIENH dataset and desikan parcellation atlas..."
## [1] "opening graphs for MRN1313 dataset and desikan parcellation atlas..."
## [1] "opening graphs for Templeton114 dataset and desikan parcellation atlas..."
## [1] "opening graphs for Templeton255 dataset and desikan parcellation atlas..."
## [1] "opening graphs for SWU4 dataset and desikan parcellation atlas..."
dwi.graphs = graphobj$graphs
dwi.datasets = graphobj$dataset
dwi.subjects = graphobj$subjects
ne = 1225
nroi <- 70
group1 <- c()  # edges in same hemisphere
group2 <- c()  # edges across hemispheres
for (i in 1:nroi) {
  for (j in 1:nroi) {
    idx <- (i - 1)*nroi + j
    if (abs(j - i) == 35) {  # across hemispheric edges
      group1 <- c(group1, idx)
    } else if (i != j) {  # ignore diagonal
      group2 <- c(group2, idx)
    }
  }
}
Es <- list(group1, group2)
dwi.models <- sapply(1:dim(dwi.graphs)[1], function(i) {
                  gs.siem.fit(thresh_matrix(dwi.graphs[i,,], 0), Es, alt='greater')
                }, simplify = FALSE)

dwi.homo.phat <- sapply(dwi.models, function(model) model$pr[1])
dwi.hetero.phat <- sapply(dwi.models, function(model)model$pr[2])
dwi.delta.phat <- sapply(dwi.models, function(model) model$dpr[1,2])
dwi.delta.var <- sapply(dwi.models, function(model)model$dvar[1,2])
dwi.phat.mu = mean(dwi.delta.phat)
dwi.phat.var = model.var(mean(dwi.homo.phat), ne) + model.var(mean(dwi.hetero.phat), ne)

Functional

nroi <- 70
fmri.dsets = c('BNU1', 'BNU2', 'BNU3', 'HNU1', 'IBATRT', 'IPCAS1', 'IPCAS2', 'IPCAS5', 'IPCAS6', 'IPCAS8', 'MRN1', 'NYU1', 'SWU1', 'SWU2', 'SWU3', 'SWU4', 'UWM', 'XHCUMS')
fmri.atlas = 'desikan-2mm'
fmri.basepath = '/data/all_mr/fmri/ranked/edgelists/'

graphobj = fmriu.io.collection.open_graphs(basepath = fmri.basepath, atlases = fmri.atlas, datasets=fmri.dsets, fmt='edgelist', rtype = 'array')
## [1] "opening graphs for BNU1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for BNU2 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for BNU3 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for HNU1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IBATRT dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS2 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS5 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS6 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS8 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for MRN1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for NYU1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for SWU1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for SWU2 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for SWU3 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for SWU4 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for UWM dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for XHCUMS dataset and desikan-2mm parcellation atlas..."
fmri.graphs = graphobj$graphs
fmri.datasets = graphobj$dataset
fmri.subjects = graphobj$subjects
ne = 1225
nroi <- 70
group1 <- c()  # edges in same hemisphere
group2 <- c()  # edges across hemispheres
for (i in 1:nroi) {
  for (j in 1:nroi) {
    idx <- (i - 1)*nroi + j
    if (abs(j - i) == 35) {  # across hemispheric edges
      group1 <- c(group1, idx)
    } else if (i != j) {  # ignore diagonal
      group2 <- c(group2, idx)
    }
  }
}
Es <- list(group1, group2)
fmri.models <- sapply(1:dim(fmri.graphs)[1], function(i) {
                  gs.siem.fit(thresh_matrix(fmri.graphs[i,,], 0.5), Es, alt='greater')
                }, simplify = FALSE)

fmri.homo.phat <- sapply(fmri.models, function(model) model$pr[1])
fmri.hetero.phat <- sapply(fmri.models, function(model)model$pr[2])
fmri.delta.phat <- sapply(fmri.models, function(model) model$dpr[1,2])
fmri.delta.var <- sapply(fmri.models, function(model)model$dvar[1,2])
fmri.phat.mu = mean(fmri.delta.phat)
fmri.phat.var = model.var(mean(fmri.homo.phat), ne) + model.var(mean(fmri.hetero.phat), ne)

\(n\) sample test

Here, we take each functional and diffusion connectome individually, and compute the parameters of our block model for each connectome. The question we seek to first answer is, given a large number of observations of \(\hat{\delta}\), can we detect a difference in homotopic vs heterotopic connectivity in the diffusion connectomes compared to the functional connectomes?

We might want to visualize the distribution of \(\delta = \hat{p}_{homo} - \hat{p}_{hetero}\) under the analytical model and compare to our empirical model for functional and diffusion:

ne = 1225
# density estimates of the two populations of delta
dwi.distr.emp.mod = density(as.numeric(dwi.delta.phat))
fmri.distr.emp.mod = density(as.numeric(fmri.delta.phat))

# variances sum for analytical computation
dwi.distr.ana = dnorm(dwi.distr.emp.mod$x, mean=dwi.phat.mu, sd=sqrt(dwi.phat.var))
fmri.distr.ana = dnorm(fmri.distr.emp.mod$x, mean=fmri.phat.mu, sd=sqrt(fmri.phat.var))

n_diff = length(dwi.distr.emp.mod$x)
dwi.dat = data.frame(x = c(dwi.distr.emp.mod$x, dwi.distr.emp.mod$x), y = c(dwi.distr.emp.mod$y, dwi.distr.ana),
                      distribution=c(rep("empirical", n_diff), rep("analytical", n_diff)))
dwi.dat$distribution = factor(dwi.dat$distribution)

ggplot(dat=dwi.dat, aes(x=x, y=y, ymax=y, fill=distribution, color=distribution, group=distribution)) +
  geom_ribbon(ymin=0, alpha=0.5) +
  ylab('Density') +
  xlab(TeX('$\\delta_D$')) +
  ggtitle(TeX('Distribution of $\\delta_D = \\hat{p}_{homo} - \\hat{p}_{hetero}$, DWI')) +
  theme(panel.background = element_rect(fill = '#ffffff'))

n_diff = length(dwi.distr.emp.mod$x)
fmri.dat = data.frame(x = c(fmri.distr.emp.mod$x, fmri.distr.emp.mod$x), y = c(fmri.distr.emp.mod$y, fmri.distr.ana),
                      distribution=c(rep("empirical", n_diff), rep("analytical", n_diff)))
fmri.dat$distribution = factor(fmri.dat$distribution)

ggplot(dat=fmri.dat, aes(x=x, y=y, ymax=y, fill=distribution, color=distribution, group=distribution)) +
  geom_ribbon(ymin=0, alpha=0.5) +
  ylab('Density') +
  xlab(TeX('$\\delta_F$')) +
  ggtitle(TeX('Distribution of $\\delta_F = \\hat{p}_{homo} - \\hat{p}_{hetero}$, fMRI')) +
  theme(panel.background = element_rect(fill = '#ffffff'))

Next, we look at the distributions, empirically and analytically, to see if there is a visually perceptible difference in the distribution of the populations \(\delta_F\) and \(\delta_D\) are from:

cmp.emp = data.frame(x = c(dwi.distr.emp.mod$x, fmri.distr.emp.mod$x), y = c(dwi.distr.emp.mod$y, fmri.distr.emp.mod$y),
                     distribution=c(rep("DWI", n_diff), rep("fMRI", n_diff)))

ggplot(dat=cmp.emp, aes(x=x, y=y, ymax=y, fill=distribution, color=distribution, group=distribution)) +
  geom_ribbon(ymin=0, alpha=0.5) +
  ylab('Density') +
  xlab(TeX('$\\delta$')) +
  ggtitle(TeX('Distribution of $\\delta$, Empirical Comparison')) +
  theme(panel.background = element_rect(fill = '#ffffff'))

cmp.ana = data.frame(x = c(dwi.distr.emp.mod$x, fmri.distr.emp.mod$x), y = c(dwi.distr.ana, fmri.distr.ana),
                     distribution=c(rep("DWI", n_diff), rep("fMRI", n_diff)))

ggplot(dat=cmp.ana, aes(x=x, y=y, ymax=y, fill=distribution, color=distribution, group=distribution)) +
  geom_ribbon(ymin=0, alpha=0.5) +
  ylab('Density') +
  xlab(TeX('$\\delta$')) +
  ggtitle(TeX('Distribution of $\\delta$, Analytical Comparison')) +
  theme(panel.background = element_rect(fill = '#ffffff'))

both the empirical and analytical results show clear separation of \(\delta_D\) and \(\delta_F\). Performing a \(t-\)test of this observation, we see:

t.test(fmri.delta.phat, dwi.delta.phat, alternative="greater", var.equal=FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  fmri.delta.phat and dwi.delta.phat
## t = 167.8, df = 4630.1, p-value < 2.2e-16
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  0.2463029       Inf
## sample estimates:
## mean of x mean of y 
## 0.4879711 0.2392294

which shows that with \(p < 2.2 \times 10^{-16}\), the difference in connectivity homotopically vs. heterotopically in functional connectomes exceeds the difference in connectivity homotopically vs heterotopically in diffusion connectomes.

\(2\)-sample test

Below, we compute a \(p-\)value given 2 graphs, 1 functional and 1 diffusion, from the same subject, where we take the cartesian product of the functional connectomes with the diffusion connectomes for each subject. For example, if we had \(2\) functional connectomes and \(2\) diffusion connectomes, we would end up with \(4\) \(p-\)values, with the first functional connectome \(\delta\) compared with the first diffusion connectome \(\delta\) (\(11\) comparison), and so on for \(12\) (first functional connectome with second diffusion connectome), \(21\), \(22\). The question we seek to first answer is, given a large number of observations of \(\hat{\delta}\), can we detect a difference in homotopic vs heterotopic connectivity in the functional connectomes compared to the diffusion connectomes?

# find the common datasets
common.dsets = fmri.dsets[which(fmri.dsets %in% dwi.dsets)]
# find the common subjects
common.sub.sid = fmri.subjects[which(fmri.subjects %in% dwi.subjects)]
common.sub.did = fmri.datasets[which(fmri.subjects %in% dwi.subjects)]

per.p <- c()
per.dataset <- c()
per.subject <- c()
fmri.failed <- c()
dwi.failed <- c()
p.failed <- c()
per.dwi.delta <- c()
per.fmri.delta <- c()
unique_subs = unique(common.sub.sid)
for (i in 1:length(unique_subs)) {
  sub = unique_subs[i]
  dset = unique(fmri.datasets[which(fmri.subjects == sub)])
  fmri.idxs = which(fmri.subjects == sub)
  dwi.idxs = which(dwi.subjects == sub)
  for (fidx in fmri.idxs) {
    for (didx in dwi.idxs) {
      pval <- gs.siem.sample.test(fmri.delta.phat[fidx], dwi.delta.phat[didx], fmri.delta.var[fidx],
                                  dwi.delta.var[didx], alt='greater', df=2)$p
      per.dwi.delta <- c(per.dwi.delta, dwi.delta.phat[didx])
      per.fmri.delta <- c(per.fmri.delta, fmri.delta.phat[fidx])
      per.p <- c(per.p, pval)
      per.dataset <- c(per.dataset, dset)
      per.subject <- c(per.subject, sub)
      if (pval > .05) {
        fmri.failed <- c(fmri.failed, fidx)
        dwi.failed <- c(dwi.failed, didx)
        p.failed <- c(p.failed, pval)
      }
    }
  }
}

p.dat <- data.frame(p = per.p, dataset=per.dataset, subject=per.subject)
p.dat$dataset <- factor(p.dat$dataset)
ggplot(data=p.dat, aes(x=dataset, y=p, color=dataset, group=dataset)) +
  geom_jitter() +
  coord_trans(y = "log10") +
  ggtitle(TeX(sprintf('Hemispheric Inter-Modality, %.2f percent have $p < .05$', 100*sum(per.p < .05)/length(per.p)))) +
  xlab('Dataset') +
  ylab('p-value') +
  theme(axis.text.x = element_text(angle=45), legend.position=NaN, panel.background = element_rect(fill = '#ffffff'))

As we can see, even at just the \(2-\)graph level, the difference in homotopic vs. heterotopic connectivity in the functional connectomes exceeds that of the diffusion connectomes and is significant in most connectomes at \(\alpha=0.05\).

We can instead visualize this as density estimates:

vline = data.frame(x=.05, type="sig")
ggplot(data=p.dat, aes(p, group=dataset, color=dataset)) +
  geom_line(stat="density", size=1.5, adjust=1.5) +
  scale_x_log10(limits=c(.005, 1)) +
  geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
  scale_color_discrete(name="Dataset") +
  scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) +
  xlab(TeX('$log(p)$')) +
  ylab("Density") +
  theme(panel.background = element_rect(fill = '#ffffff')) + 
  ggtitle(TeX(sprintf('Bilateral Inter-Modality, %.2f percent have $p < .05$', 100*sum(p.dat$p < .05)/length(p.dat$p))))

Example of subject with \(p > .05\)
sort_p <- sort(p.failed, index.return=TRUE, decreasing = TRUE)
ix = sort_p$ix[1]
dwfailed = thresh_matrix(dwi.graphs[dwi.failed[ix],,], thresh=0)
fmfailed = thresh_matrix(fmri.graphs[fmri.failed[ix],,])
fmriu.plot.plot_graph(dwfailed, include_diag = TRUE,
                      title = sprintf("DWI subject %s, delta=%.3f", as.character(dwi.subjects[dwi.failed[ix]]), dwi.delta.phat[dwi.failed[ix]]),
                      legend.name = "connection")

fmriu.plot.plot_graph(fmfailed, include_diag = TRUE,
                      title = sprintf("fMRI subject %s, delta=%.3f", as.character(fmri.subjects[fmri.failed[ix]]), fmri.delta.phat[fmri.failed[ix]]),
                      legend.name = "connection")

print(sort_p$x[1])
## [1] 0.9523378
sort_p <- sort(p.failed, index.return=TRUE, decreasing = TRUE)
ix = sort_p$ix[2]
dwfailed = thresh_matrix(dwi.graphs[dwi.failed[ix],,], thresh=0)
fmfailed = thresh_matrix(fmri.graphs[fmri.failed[ix],,])
fmriu.plot.plot_graph(dwfailed, include_diag = TRUE,
                      title = sprintf("DWI subject %s, delta=%.3f", as.character(dwi.subjects[dwi.failed[ix]]), dwi.delta.phat[dwi.failed[ix]]),
                      legend.name = "connection")

fmriu.plot.plot_graph(fmfailed, include_diag = TRUE,
                      title = sprintf("fMRI subject %s, delta=%.3f", as.character(fmri.subjects[fmri.failed[ix]]), fmri.delta.phat[fmri.failed[ix]]),
                      legend.name = "connection")

print(sort_p$x[2])
## [1] 0.9503484

Below, we visualize the difference \(\delta_D - \delta_F\) between each pair of connectomes, and investigate the result in the \(p-\)value:

batch.dat <- data.frame(diff= per.dwi.delta - per.fmri.delta, p = per.p, dataset=per.dataset)
batch.dat$dataset <- factor(batch.dat$dataset)
ggplot(batch.dat, aes(x=diff, y=p, color=dataset, group=dataset, shape=dataset)) +
  geom_point() +
  xlab(TeX('$\\delta_D - \\delta_F$')) +
  ylab(TeX('$p-$value')) +
  ggtitle(TeX('Examining impact on $p-$value from difference in connectivity')) +
  theme(panel.background = element_rect(fill = '#ffffff'))

Megamean

Here, we again perform a test on 2 graphs, except here the graphs used are the average functional and diffusion connectomes (the megameans). We feed this into a simple t-test with the appropriate assumptions (unequal variance, goal is to test for homotopic connectivity exceeding heterotopic connectivity):

fmri.agg.mod <- gs.siem.fit(fmri_thresh, Es)
dwi.agg.mod <- gs.siem.fit(dwi_thresh, Es)

gs.siem.sample.test(fmri.agg.mod$dpr[1,2], dwi.agg.mod$dpr[1,2], fmri.agg.mod$dvar[1,2], dwi.agg.mod$dvar[1,2], df = 2)$p
## [1] 0.01336706

As we can see, with just 2 megamean graphs, we can identify a significant difference in the difference in connectivity homotopically vs bi-laterally in functional vs. diffusion connectomes.