require(fmriutils)
require(graphstats)
require(mgc)
require(ggplot2)
require(latex2exp)
require(igraph)
source('../../R/siem.R')

# accepts a matrix and thresholds/binarizes it
thresh_matrix = function(matrix, thresh=0.5) {
  thr = quantile(matrix, thresh)
  return(ifelse(matrix > thr, 1, 0))
}

# accepts a [n x n] adjacency matrix and computes the probabilities associated with an SBM
# where the vertices are grouped such that union_i(groups) = V(A) and
# intersection(group_i, group_j) = NULL for all i != j
block_data = function(matrix, groups) {
  # matrix is adwi_thresh n x n array
  # groups is a grouping of the vertices in the matrix as a list
  blocks = array(NaN, dim=c(2,2))
  for (i in 1:length(groups)) {
    for (j in 1:length(groups)) {
      blocks[i, j] = mean(matrix[groups[[i]], groups[[j]]])
    }
  }
  return(blocks)
}

Model

For details on the model, please see SIEM model.

Real Data Experiments

In this notebook, we determine whether there exists a significant difference in connectivity ipsi-laterally (same hemisphere) vs contra-laterally (opposite hemisphere) connectivity within a particular modality. We consider $\delta_{x} = p_{ipsi} - p_{contra}$ 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 ipsi-laterally exceeds the connectivity contra-laterally. We will do this in 2 ways:

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 ipsilateral quadrants; 2 contralateral 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')
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 ((i <= 35 & j <= 35) | (i > 35 & j > 35)) {
      group1 <- c(group1, idx)
    } else {
      group2 <- c(group2, idx)
    }
  }
}
Es <- list(group1, group2)
dwi.models <- sapply(1:dim(dwi.graphs)[1], function(i) {
                  list(model=gs.siem.fit(thresh_matrix(dwi.graphs[i,,], 0), Es, alt='greater'), Z=dwi.datasets[i])
                }, simplify = FALSE)

dwi.ips.phat <- sapply(dwi.models, function(model) model$model$pr[1])
dwi.contr.phat <- sapply(dwi.models, function(model)model$model$pr[2])
dwi.ips.phat <- sapply(dwi.models, function(model) model$model$pr[1])
dwi.contr.phat <- sapply(dwi.models, function(model)model$model$pr[2])
dwi.delta.phat <- sapply(dwi.models, function(model) model$model$dpr[1,2])
dwi.delta.var <- sapply(dwi.models, function(model)model$model$dvar[1,2])
dwi.phat.mu = mean(dwi.delta.phat)
dwi.phat.var = model.var(mean(dwi.ips.phat), ne) + model.var(mean(dwi.contr.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')
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 ((i <= 35 & j <= 35) | (i > 35 & j > 35)) {
      group1 <- c(group1, idx)
    } else {
      group2 <- c(group2, idx)
    }
  }
}
Es <- list(group1, group2)
fmri.models <- sapply(1:dim(fmri.graphs)[1], function(i) {
                  list(model=gs.siem.fit(thresh_matrix(fmri.graphs[i,,], 0.5), Es, alt='greater'), Z=dwi.datasets[i])
                }, simplify = FALSE)

fmri.ips.phat <- sapply(fmri.models, function(model) model$model$pr[1])
fmri.contr.phat <- sapply(fmri.models, function(model)model$model$pr[2])
fmri.delta.phat <- sapply(fmri.models, function(model) model$model$dpr[1,2])
fmri.delta.var <- sapply(fmri.models, function(model)model$model$dvar[1,2])
fmri.phat.mu = mean(fmri.delta.phat)
fmri.phat.var = model.var(mean(fmri.ips.phat), ne) + model.var(mean(fmri.contr.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 ipsi-lateral vs contra-lateral connectivity in the diffusion connectomes compared to the functional connectomes?

We might want to visualize the distribution of $\delta = |\hat{p}{contra} - \hat{p}{ipsi}|$ 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}_{ipsi} - \\hat{p}_{contr}|$, 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}_{ipsi} - \\hat{p}_{contr}|$, 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(dwi.delta.phat, fmri.delta.phat, alternative="greater", var.equal=FALSE)

which shows that with $p < 2.2 \times 10^{-16}$, the difference in connectivity ipsi-laterally vs. contra-laterally in diffusion connectomes exceeds the difference in connectivity ipsi-laterally vs contra-laterally in functional 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 ipsi-lateral vs contra-lateral 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(dwi.delta.phat[didx], fmri.delta.phat[fidx], dwi.delta.var[didx],
                                  fmri.delta.var[fidx], 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 ipsi-lateral vs. contra-lateral connectivity in the diffusion connectomes exceeds that of the functional 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('Hemispheric 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])
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])

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'))

Aggregated

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 ipsilateral connectivity exceeding contralateral connectivity):

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

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

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



neurodata/graphstats documentation built on May 14, 2019, 5:19 p.m.