In neurodata/graphstats: Graph Statistics
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 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')
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')
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)
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])
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'))
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
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.
neurodata/graphstats documentation built on May 14, 2019, 5:19 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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 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') 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') 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)
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])
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'))
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
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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.