In neurodata/graphstats: Graph Statistics
require(fmriutils)
require(graphstats)
require(ggplot2)
require(latex2exp)
require(igraph)
require(scales)
require(gridExtra)
require(grid)
require(data.table)
is.defined = function(x)!is.null(x)
# accepts a matrix and thresholds/binarizes it
thresh_matrix = function(matrix, thresh=0.5) {
thr = quantile(matrix, thresh)
return(ifelse(matrix > thr, 1, 0))
}
compute_params <- function(graphs, datasets, subjects, edges, thresh=-1) {
if (thresh >= 0) {
graphs <- lapply(1:dim(graphs)[1], function(i) thresh_matrix(graphs[i,,], thresh))
}
models <- suppressWarnings(lapply(graphs, function(graph) gs.siem.fit(graph, edges)))
incl <- which(sapply(models, is.defined))
models <- models[incl]
dsets <- datasets[incl]
subs <- subjects[incl]
e1.phat <- sapply(models, function(model) model$pr[1])
e2.phat <- sapply(models, function(model)model$pr[2])
delta <- sapply(models, function(model) model$dpr[1,2])
delta.var <- sapply(models, function(model)model$dvar[1,2])
phat.mu <- mean(delta)
phat.var <- graphstats:::model.var(mean(e1.phat), ne) + graphstats:::model.var(mean(e2.phat), ne)
return(list(models=models, datasets=datasets, subjects=subs, e1.phat=e1.phat, e2.phat=e2.phat,
delta=delta, phat.mu=phat.mu, phat.var=phat.var))
}
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)
}
}
}
bilateral.edges <- list(group1, group2)
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)
}
}
}
hemisphere.edges <- list(group1, group2)
vline = data.frame(x=.05, type="sig")
model.between <- function(fit) {
model.diff <- lapply(1:length(fit$models), function(i) {
diff = fit$models[[i]]$pr[1] - fit$models[[i]]$pr[2]
pval <- fit$models[[i]]$pv[1, 2]
data.frame(difference = diff, dataset=fit$datasets[i], pval=pval, dpr=fit$models[[i]]$dpr[1,2],
dvar=fit$models[[i]]$dvar[1,2], subject=fit$subjects[i])
})
return(do.call(rbind, model.diff))
}
model.across <- function(within.first, within.second) {
first.subs <- within.first$subject; second.subs <- within.second$subject
subs <- first.subs[first.subs %in% second.subs]
model.acr <- lapply(1:length(subs), function(i) {
sec <- within.second[as.character(second.subs) == subs[i],]
fir <- within.first[as.character(first.subs) == subs[i],]
res.out <- lapply(1:nrow(fir), function(j) {
sub.first <- fir[j,]
res.in <- lapply(1:nrow(sec), function(k) {
sub.second <- sec[k,]
pval <- gs.siem.sample.test(fir$dpr, sec$dpr, fir$dvar, sec$dvar, alt='greater', df=2)$p
data.frame(first = sub.first$difference, second=sub.second$difference,
difference = sub.first$difference - sub.second$difference, pval=pval,
dataset = as.character(sub.first$dataset), subject=as.character(subs[i]))
})
do.call(rbind, res.in)
})
do.call(rbind, res.out)
})
return(do.call(rbind, model.acr))
}
g_legend <- function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)
}
colors <- readRDS('./data/colors.rds')
svals <- colors$svals; cols <- colors$cols
svals['BNU'] <- svals['BNU1'] + svals['BNU3']
svals['Templeton'] <- svals['Templeton114'] + svals['Templeton255']
svals['SWU'] <- svals['SWU1'] + svals['SWU2'] + svals['SWU3'] + svals['SWU4']
svals['IPCAS'] <- svals['IPCAS1'] + svals['IPCAS2'] + svals['IPCAS5'] + svals['IPCAS6'] + svals['IPCAS8']
cols['SWU'] <- '#5f5a0c'; cols['IPCAS'] <- '#ff0000'
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
dwi.hem.fit <- compute_params(dwi.graphs, dwi.datasets, dwi.subjects, hemisphere.edges, thresh=0)
dwi.bil.fit <- compute_params(dwi.graphs, dwi.datasets, dwi.subjects, bilateral.edges, thresh=0)
Functional
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
fmri.hem.fit <- compute_params(fmri.graphs, fmri.datasets, fmri.subjects, hemisphere.edges, thresh=0.5)
fmri.bil.fit <- compute_params(fmri.graphs, fmri.datasets, fmri.subjects, bilateral.edges, thresh=0.5)
Hemisphieric
Within Modality
dMRI
dwi.hem.diff <- model.between(dwi.hem.fit)
dwi.hem.diff$modality = 'dMRI'
fMRI
fmri.hem.diff <- model.between(fmri.hem.fit)
fmri.hem.diff$modality = 'fMRI'
Plot
vline = data.frame(x=0.0, type="sig")
hem.diff <- rbind(dwi.hem.diff, fmri.hem.diff)
hem.diff$modality <- factor(hem.diff$modality, levels=unique(hem.diff$modality))
hem.diff$grouping = paste(hem.diff$dataset, hem.diff$modality) # for the datasets that are shared
within.hem <- ggplot(data=hem.diff, aes(difference, group=grouping, color=modality)) +
geom_line(aes(y=..scaled..), stat="density", size=1, adjust=1.5) +
xlab(TeX('$\\delta = \\hat{p} - \\hat{q}$')) +
geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
scale_linetype_manual(values=c("dotted"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\delta_D = \\delta_F$"), TeX)) +
xlim(-1, 1) +
ylab("Normalized Density") +
theme_bw() +
theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) +
labs(color="Modality") +
guides(linetype=FALSE) +
ggtitle("Hemispheric Intra-Modality")
CDFs
vline = data.frame(x=0.05, type="sig")
hem.cdf <- ggplot(data=hem.diff, aes(x=pval, color=modality)) +
stat_ecdf(size=1, geom="line") +
scale_x_log10(limits=c(.001, 1), breaks=c(.001, .01, .1, 1), labels=c(.001, .01, .1, 1)) +
geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) +
xlab(TeX('$p$-value')) +
ylab("Cumulative Distribution") +
theme_bw() +
theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) +
labs(color="Modality") +
guides(color=FALSE) +
ggtitle("Cumulative Distribution of P-values")
Within Modality
Plot
Bilateral
Within Modality
dMRI
dwi.bil.diff <- model.between(dwi.bil.fit)
dwi.bil.diff$modality = 'dMRI'
fMRI
fmri.bil.diff <- model.between(fmri.bil.fit)
fmri.bil.diff$modality = 'fMRI'
Plot
vline = data.frame(x=0.0, type="sig")
bil.diff <- rbind(dwi.bil.diff, fmri.bil.diff)
bil.diff$modality <- factor(bil.diff$modality, levels=unique(bil.diff$modality))
bil.diff$grouping = paste(bil.diff$dataset, bil.diff$modality) # for the datasets that are shared
within.bil <- ggplot(data=bil.diff, aes(difference, group=grouping, color=modality)) +
geom_line(aes(y=..scaled..), stat="density", size=1, adjust=1.5) +
xlab(TeX('$\\delta = \\hat{p} - \\hat{q}$')) +
geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
scale_linetype_manual(values=c("dotted"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\delta_D = \\delta_F$"), TeX)) +
ylab("Normalized Density") +
theme_bw() +
theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) +
xlim(-1, 1) +
labs(color="Modality") +
guides(linetype=FALSE) +
ggtitle("Bilateral Intra-Modality")
CDFs
vline = data.frame(x=0.05, type="sig")
bil.cdf <- ggplot(data=bil.diff, aes(x=pval, color=modality)) +
stat_ecdf(size=1, geom="line") +
scale_x_log10(limits=c(.001, 1), breaks=c(.001, .01, .1, 1), labels=c(.001, .01, .1, 1)) +
geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) +
xlab(TeX('$p$-value')) +
ylab("Cumulative Distribution") +
theme_bw() +
theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) +
labs(color="Modality") +
guides(color=FALSE) +
ggtitle("Cumulative Distribution of P-values")
Between Modality
hem.across <- model.across(dwi.hem.diff, fmri.hem.diff)
bil.across <- model.across(fmri.bil.diff, dwi.bil.diff)
Plot
vline = data.frame(x=0, type="sig")
across.hem <- ggplot(data=hem.across, aes(difference, color=dataset)) +
geom_line(aes(y=..scaled..), stat="density", size=1, adjust=1.5) +
geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
scale_linetype_manual(values=c("dotted"), name="Cutoff", breaks=c("sig"), labels=lapply(c("ipsi = contra"), TeX)) +
xlim(-1, 1) +
xlab(TeX('$\\hat{\\delta}_D - \\hat{\\delta}_F')) +
ylab("Normalized Density") +
scale_color_manual(values=cols) +
theme_bw() +
theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) +
labs(color="Dataset") +
guides(linetype=FALSE) +
ggtitle("Hemispheric Across Modality")
vline = data.frame(x=0, type="sig")
across.bil <- ggplot(data=bil.across, aes(difference, color=dataset)) +
geom_line(aes(y=..scaled..), stat="density", size=1, adjust=1.5) +
geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
scale_linetype_manual(values=c("dotted"), name="Cutoff", breaks=c("sig"), labels=lapply(c("homo = hetero"), TeX)) +
xlim(-1, 1) +
xlab(TeX('$\\hat{\\delta}_F - \\hat{\\delta}_D')) +
ylab("Normalized Density") +
scale_color_manual(values=cols) +
theme_bw() +
theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) +
labs(color="Dataset") +
guides(linetype=FALSE) +
ggtitle("Bilateral Across Modality")
CDFs
vline = data.frame(x=0.05, type="sig")
across.hem.cdf <- ggplot(data=hem.across, aes(x=pval, color=dataset)) +
stat_ecdf(size=1, geom="line") +
scale_color_manual(values=cols) +
scale_x_log10(limits=c(.001, 1), breaks=c(.001, .01, .1, 1), labels=c(.001, .01, .1, 1)) +
geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) +
xlab(TeX('$p$-value')) +
ylab("Cumulative Distribution") +
theme_bw() +
theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) +
labs(color="Modality") +
guides(color=FALSE) +
ggtitle("Cumulative Distribution of P-values")
across.bil.cdf <- ggplot(data=bil.across, aes(x=pval, color=dataset)) +
stat_ecdf(size=1, geom="line") +
scale_x_log10(limits=c(.001, 1), breaks=c(.001, .01, .1, 1), labels=c(.001, .01, .1, 1)) +
geom_vline(data=vline, aes(xintercept = x, linetype=type)) +
scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) +
xlab(TeX('$p$-value')) +
scale_color_manual(values=cols) +
ylab("Cumulative Distribution") +
theme_bw() +
theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) +
labs(color="Modality") +
guides(color=FALSE) +
ggtitle("Cumulative Distribution of P-values")
saveRDS(list(hem.across=hem.across, bil.across=bil.across, hem.diff=hem.diff, bil.diff=bil.diff), 'data/figure_4.rds')
Combined
c1_leg <- g_legend(within.hem + guides(color = guide_legend(nrow = 1)))
c2_leg <- g_legend(hem.cdf + guides(linestyle = guide_legend(nrow = 1)))
c3_leg <- g_legend(across.hem + guides(color = guide_legend(nrow = 1)))
c4_leg <- g_legend(across.hem.cdf + guides(linestyle = guide_legend(nrow = 1)))
r_titls <- list(r1_titl="Hemispheric", r2_titl="Homotopic")
c_titls <- list(c1_titl="Magnitude Within-Modality", c2_titl="CDF Within-Modality", c3_titl="Magnitude Inter-Modality", c4_titl="CDF Inter-Modality")
col.titles <- lapply(c_titls, function(titl) textGrob(titl, gp=gpar(fontsize=26)))
row.titles <- lapply(r_titls, function(titl) textGrob(titl, gp=gpar(fontsize=15), rot=90))
empty.title <- textGrob("", gp=gpar(fontsize=18))
col1 <- arrangeGrob(within.hem + ggtitle("(A.i)") + xlab("ipsi - contra") + ylab("Normalized Density") + theme(legend.position=NaN),
within.bil + ggtitle("(B.i)") + xlab("homo - hetero") + ylab("") + theme(legend.position=NaN),
c1_leg, nrow=3, heights=c(.45, .45, .15), top=c_titls$c1_titl)
col2 <- arrangeGrob(hem.cdf + ggtitle("(A.ii)") + xlab(TeX("$p$-value")) + ylab("") + theme(legend.position=NaN),
bil.cdf + ggtitle("(B.ii)") + xlab("") + ylab("") + theme(legend.position=NaN),
c2_leg, nrow=3, heights=c(.45, .45, .15), top=c_titls$c2_titl)
col3 <- arrangeGrob(across.hem + ggtitle("(A.iii)") + xlab("dMRI - fMRI") + ylab("") + theme(legend.position=NaN),
across.bil + ggtitle("(B.iii)") + xlab("fMRI - dMRI") + ylab("") + theme(legend.position=NaN),
c3_leg, nrow=3, heights=c(.45, .45, .15), top=c_titls$c3_titl)
col4 <- arrangeGrob(across.hem.cdf + ggtitle("(A.iv)") + xlab(TeX("$p$-value")) + ylab("") + theme(legend.position=NaN),
across.bil.cdf + ggtitle("(B.iv)") + xlab("") + ylab("") + theme(legend.position=NaN),
c4_leg, nrow=3, heights=c(.45, .45, .15), top=c_titls$c4_titl)
col_left <- arrangeGrob(row.titles$r1_titl, row.titles$r2_titl, empty.title, nrow=3, heights=c(.45, .45, .1))
grid.arrange(col_left, col1, col2, col3, col4, widths=c(.02, .2, .2, .2, .2))
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(ggplot2) require(latex2exp) require(igraph) require(scales) require(gridExtra) require(grid) require(data.table) is.defined = function(x)!is.null(x) # accepts a matrix and thresholds/binarizes it thresh_matrix = function(matrix, thresh=0.5) { thr = quantile(matrix, thresh) return(ifelse(matrix > thr, 1, 0)) } compute_params <- function(graphs, datasets, subjects, edges, thresh=-1) { if (thresh >= 0) { graphs <- lapply(1:dim(graphs)[1], function(i) thresh_matrix(graphs[i,,], thresh)) } models <- suppressWarnings(lapply(graphs, function(graph) gs.siem.fit(graph, edges))) incl <- which(sapply(models, is.defined)) models <- models[incl] dsets <- datasets[incl] subs <- subjects[incl] e1.phat <- sapply(models, function(model) model$pr[1]) e2.phat <- sapply(models, function(model)model$pr[2]) delta <- sapply(models, function(model) model$dpr[1,2]) delta.var <- sapply(models, function(model)model$dvar[1,2]) phat.mu <- mean(delta) phat.var <- graphstats:::model.var(mean(e1.phat), ne) + graphstats:::model.var(mean(e2.phat), ne) return(list(models=models, datasets=datasets, subjects=subs, e1.phat=e1.phat, e2.phat=e2.phat, delta=delta, phat.mu=phat.mu, phat.var=phat.var)) } 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) } } } bilateral.edges <- list(group1, group2) 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) } } } hemisphere.edges <- list(group1, group2) vline = data.frame(x=.05, type="sig") model.between <- function(fit) { model.diff <- lapply(1:length(fit$models), function(i) { diff = fit$models[[i]]$pr[1] - fit$models[[i]]$pr[2] pval <- fit$models[[i]]$pv[1, 2] data.frame(difference = diff, dataset=fit$datasets[i], pval=pval, dpr=fit$models[[i]]$dpr[1,2], dvar=fit$models[[i]]$dvar[1,2], subject=fit$subjects[i]) }) return(do.call(rbind, model.diff)) } model.across <- function(within.first, within.second) { first.subs <- within.first$subject; second.subs <- within.second$subject subs <- first.subs[first.subs %in% second.subs] model.acr <- lapply(1:length(subs), function(i) { sec <- within.second[as.character(second.subs) == subs[i],] fir <- within.first[as.character(first.subs) == subs[i],] res.out <- lapply(1:nrow(fir), function(j) { sub.first <- fir[j,] res.in <- lapply(1:nrow(sec), function(k) { sub.second <- sec[k,] pval <- gs.siem.sample.test(fir$dpr, sec$dpr, fir$dvar, sec$dvar, alt='greater', df=2)$p data.frame(first = sub.first$difference, second=sub.second$difference, difference = sub.first$difference - sub.second$difference, pval=pval, dataset = as.character(sub.first$dataset), subject=as.character(subs[i])) }) do.call(rbind, res.in) }) do.call(rbind, res.out) }) return(do.call(rbind, model.acr)) } g_legend <- function(a.gplot){ tmp <- ggplot_gtable(ggplot_build(a.gplot)) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend) }
colors <- readRDS('./data/colors.rds') svals <- colors$svals; cols <- colors$cols svals['BNU'] <- svals['BNU1'] + svals['BNU3'] svals['Templeton'] <- svals['Templeton114'] + svals['Templeton255'] svals['SWU'] <- svals['SWU1'] + svals['SWU2'] + svals['SWU3'] + svals['SWU4'] svals['IPCAS'] <- svals['IPCAS1'] + svals['IPCAS2'] + svals['IPCAS5'] + svals['IPCAS6'] + svals['IPCAS8'] cols['SWU'] <- '#5f5a0c'; cols['IPCAS'] <- '#ff0000'
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 dwi.hem.fit <- compute_params(dwi.graphs, dwi.datasets, dwi.subjects, hemisphere.edges, thresh=0) dwi.bil.fit <- compute_params(dwi.graphs, dwi.datasets, dwi.subjects, bilateral.edges, thresh=0)
Functional
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 fmri.hem.fit <- compute_params(fmri.graphs, fmri.datasets, fmri.subjects, hemisphere.edges, thresh=0.5) fmri.bil.fit <- compute_params(fmri.graphs, fmri.datasets, fmri.subjects, bilateral.edges, thresh=0.5)
Hemisphieric
Within Modality
dMRI
dwi.hem.diff <- model.between(dwi.hem.fit) dwi.hem.diff$modality = 'dMRI'
fMRI
fmri.hem.diff <- model.between(fmri.hem.fit) fmri.hem.diff$modality = 'fMRI'
Plot
vline = data.frame(x=0.0, type="sig") hem.diff <- rbind(dwi.hem.diff, fmri.hem.diff) hem.diff$modality <- factor(hem.diff$modality, levels=unique(hem.diff$modality)) hem.diff$grouping = paste(hem.diff$dataset, hem.diff$modality) # for the datasets that are shared within.hem <- ggplot(data=hem.diff, aes(difference, group=grouping, color=modality)) + geom_line(aes(y=..scaled..), stat="density", size=1, adjust=1.5) + xlab(TeX('$\\delta = \\hat{p} - \\hat{q}$')) + geom_vline(data=vline, aes(xintercept = x, linetype=type)) + scale_linetype_manual(values=c("dotted"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\delta_D = \\delta_F$"), TeX)) + xlim(-1, 1) + ylab("Normalized Density") + theme_bw() + theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) + labs(color="Modality") + guides(linetype=FALSE) + ggtitle("Hemispheric Intra-Modality")
CDFs
vline = data.frame(x=0.05, type="sig") hem.cdf <- ggplot(data=hem.diff, aes(x=pval, color=modality)) + stat_ecdf(size=1, geom="line") + scale_x_log10(limits=c(.001, 1), breaks=c(.001, .01, .1, 1), labels=c(.001, .01, .1, 1)) + geom_vline(data=vline, aes(xintercept = x, linetype=type)) + scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) + xlab(TeX('$p$-value')) + ylab("Cumulative Distribution") + theme_bw() + theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) + labs(color="Modality") + guides(color=FALSE) + ggtitle("Cumulative Distribution of P-values")
Within Modality
Plot
Bilateral
Within Modality
dMRI
dwi.bil.diff <- model.between(dwi.bil.fit) dwi.bil.diff$modality = 'dMRI'
fMRI
fmri.bil.diff <- model.between(fmri.bil.fit) fmri.bil.diff$modality = 'fMRI'
Plot
vline = data.frame(x=0.0, type="sig") bil.diff <- rbind(dwi.bil.diff, fmri.bil.diff) bil.diff$modality <- factor(bil.diff$modality, levels=unique(bil.diff$modality)) bil.diff$grouping = paste(bil.diff$dataset, bil.diff$modality) # for the datasets that are shared within.bil <- ggplot(data=bil.diff, aes(difference, group=grouping, color=modality)) + geom_line(aes(y=..scaled..), stat="density", size=1, adjust=1.5) + xlab(TeX('$\\delta = \\hat{p} - \\hat{q}$')) + geom_vline(data=vline, aes(xintercept = x, linetype=type)) + scale_linetype_manual(values=c("dotted"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\delta_D = \\delta_F$"), TeX)) + ylab("Normalized Density") + theme_bw() + theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) + xlim(-1, 1) + labs(color="Modality") + guides(linetype=FALSE) + ggtitle("Bilateral Intra-Modality")
CDFs
vline = data.frame(x=0.05, type="sig") bil.cdf <- ggplot(data=bil.diff, aes(x=pval, color=modality)) + stat_ecdf(size=1, geom="line") + scale_x_log10(limits=c(.001, 1), breaks=c(.001, .01, .1, 1), labels=c(.001, .01, .1, 1)) + geom_vline(data=vline, aes(xintercept = x, linetype=type)) + scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) + xlab(TeX('$p$-value')) + ylab("Cumulative Distribution") + theme_bw() + theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) + labs(color="Modality") + guides(color=FALSE) + ggtitle("Cumulative Distribution of P-values")
Between Modality
hem.across <- model.across(dwi.hem.diff, fmri.hem.diff) bil.across <- model.across(fmri.bil.diff, dwi.bil.diff)
Plot
vline = data.frame(x=0, type="sig") across.hem <- ggplot(data=hem.across, aes(difference, color=dataset)) + geom_line(aes(y=..scaled..), stat="density", size=1, adjust=1.5) + geom_vline(data=vline, aes(xintercept = x, linetype=type)) + scale_linetype_manual(values=c("dotted"), name="Cutoff", breaks=c("sig"), labels=lapply(c("ipsi = contra"), TeX)) + xlim(-1, 1) + xlab(TeX('$\\hat{\\delta}_D - \\hat{\\delta}_F')) + ylab("Normalized Density") + scale_color_manual(values=cols) + theme_bw() + theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) + labs(color="Dataset") + guides(linetype=FALSE) + ggtitle("Hemispheric Across Modality") vline = data.frame(x=0, type="sig") across.bil <- ggplot(data=bil.across, aes(difference, color=dataset)) + geom_line(aes(y=..scaled..), stat="density", size=1, adjust=1.5) + geom_vline(data=vline, aes(xintercept = x, linetype=type)) + scale_linetype_manual(values=c("dotted"), name="Cutoff", breaks=c("sig"), labels=lapply(c("homo = hetero"), TeX)) + xlim(-1, 1) + xlab(TeX('$\\hat{\\delta}_F - \\hat{\\delta}_D')) + ylab("Normalized Density") + scale_color_manual(values=cols) + theme_bw() + theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) + labs(color="Dataset") + guides(linetype=FALSE) + ggtitle("Bilateral Across Modality")
CDFs
vline = data.frame(x=0.05, type="sig") across.hem.cdf <- ggplot(data=hem.across, aes(x=pval, color=dataset)) + stat_ecdf(size=1, geom="line") + scale_color_manual(values=cols) + scale_x_log10(limits=c(.001, 1), breaks=c(.001, .01, .1, 1), labels=c(.001, .01, .1, 1)) + geom_vline(data=vline, aes(xintercept = x, linetype=type)) + scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) + xlab(TeX('$p$-value')) + ylab("Cumulative Distribution") + theme_bw() + theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) + labs(color="Modality") + guides(color=FALSE) + ggtitle("Cumulative Distribution of P-values") across.bil.cdf <- ggplot(data=bil.across, aes(x=pval, color=dataset)) + stat_ecdf(size=1, geom="line") + scale_x_log10(limits=c(.001, 1), breaks=c(.001, .01, .1, 1), labels=c(.001, .01, .1, 1)) + geom_vline(data=vline, aes(xintercept = x, linetype=type)) + scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05$"), TeX)) + xlab(TeX('$p$-value')) + scale_color_manual(values=cols) + ylab("Cumulative Distribution") + theme_bw() + theme(plot.margin = unit(c(.1,.1,.1,.1), "cm")) + labs(color="Modality") + guides(color=FALSE) + ggtitle("Cumulative Distribution of P-values")
saveRDS(list(hem.across=hem.across, bil.across=bil.across, hem.diff=hem.diff, bil.diff=bil.diff), 'data/figure_4.rds')
Combined
c1_leg <- g_legend(within.hem + guides(color = guide_legend(nrow = 1))) c2_leg <- g_legend(hem.cdf + guides(linestyle = guide_legend(nrow = 1))) c3_leg <- g_legend(across.hem + guides(color = guide_legend(nrow = 1))) c4_leg <- g_legend(across.hem.cdf + guides(linestyle = guide_legend(nrow = 1))) r_titls <- list(r1_titl="Hemispheric", r2_titl="Homotopic") c_titls <- list(c1_titl="Magnitude Within-Modality", c2_titl="CDF Within-Modality", c3_titl="Magnitude Inter-Modality", c4_titl="CDF Inter-Modality") col.titles <- lapply(c_titls, function(titl) textGrob(titl, gp=gpar(fontsize=26))) row.titles <- lapply(r_titls, function(titl) textGrob(titl, gp=gpar(fontsize=15), rot=90)) empty.title <- textGrob("", gp=gpar(fontsize=18)) col1 <- arrangeGrob(within.hem + ggtitle("(A.i)") + xlab("ipsi - contra") + ylab("Normalized Density") + theme(legend.position=NaN), within.bil + ggtitle("(B.i)") + xlab("homo - hetero") + ylab("") + theme(legend.position=NaN), c1_leg, nrow=3, heights=c(.45, .45, .15), top=c_titls$c1_titl) col2 <- arrangeGrob(hem.cdf + ggtitle("(A.ii)") + xlab(TeX("$p$-value")) + ylab("") + theme(legend.position=NaN), bil.cdf + ggtitle("(B.ii)") + xlab("") + ylab("") + theme(legend.position=NaN), c2_leg, nrow=3, heights=c(.45, .45, .15), top=c_titls$c2_titl) col3 <- arrangeGrob(across.hem + ggtitle("(A.iii)") + xlab("dMRI - fMRI") + ylab("") + theme(legend.position=NaN), across.bil + ggtitle("(B.iii)") + xlab("fMRI - dMRI") + ylab("") + theme(legend.position=NaN), c3_leg, nrow=3, heights=c(.45, .45, .15), top=c_titls$c3_titl) col4 <- arrangeGrob(across.hem.cdf + ggtitle("(A.iv)") + xlab(TeX("$p$-value")) + ylab("") + theme(legend.position=NaN), across.bil.cdf + ggtitle("(B.iv)") + xlab("") + ylab("") + theme(legend.position=NaN), c4_leg, nrow=3, heights=c(.45, .45, .15), top=c_titls$c4_titl) col_left <- arrangeGrob(row.titles$r1_titl, row.titles$r2_titl, empty.title, nrow=3, heights=c(.45, .45, .1)) grid.arrange(col_left, col1, col2, col3, col4, widths=c(.02, .2, .2, .2, .2))
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.