In neurodata/graphstats: Graph Statistics

require(fmriutils)
require(graphstats)
require(ggplot2)
require(latex2exp)
require(igraph)
require(stringr)
require(gridExtra)
require(scales)
require(data.table)
require(grid)
require(ggbeeswarm)
require(ggridges)
source('../../R/siem.R')


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'

Aggregate plots

panel_plots <- list()
files <- c("params", "pair", "case1", "case2", "case3", "case4", "case5", "case6")
title <- c("Exploratory", "(%s.%s) Between Study Differences", "Sessions", "Individuals", "Sexes", "Sites",
           "Demographics", "All Studies")
modal_plots <- list(dMRI=list(), fMRI=list())
counters <- list(dMRI=1, fMRI=1)
case_dat <- list()
nroi <- 70
ne <- nroi^2
minpos <- c((ne/2 + 1)/2, (nroi + 1)/2)
maxpos <- c((ne + ne/2 + 1)/2, (ne + ne - nroi + 1)/2)
minp <- c(0.6, 0.5)
maxp <- c(0.8, 0.7)
cutoff = 0.05
letters=c("A", "B")
temps = c("Templeton114", "Templeton255", "Templeton")
temprep = c("Temp114", "Temp255", "Temp")
xl <- c("Ipsilateral", "Homotopic")
names(svals)[names(svals) %in% temps] = temprep[match(names(svals)[names(svals) %in% temps], temps)]
names(cols)[names(cols) %in% temps] = temprep[match(names(cols)[names(cols) %in% temps], temps)]
roman=c("i","ii", "iii")
#minp <- c(0, 0)
proportion <- data.frame(case=c(), prop=c())
for (i in 1:length(names(modal_plots))) {
  modal <- names(modal_plots)[i]
  if (modal == 'dMRI') {
    path_str <- 'hem'
  } else {
    path_str <- 'bil'
  }
  top_plots <- list()
  for (pp in 1:2) {
    fpath <- paste("data/results", path_str, sprintf('%s.rds', files[pp]), sep="_")
    param = readRDS(fpath)
    if (pp == 2) {
      results <- param[[which(modal == names(modal_plots))]]$data
      results$dset1 <- as.character(results$dset1)
      results$dset2 <- as.character(results$dset2)
      results$dset1[results$dset1 %in% temps] = temprep[match(results$dset1[results$dset1 %in% temps], temps)]
      results$dset2[results$dset2 %in% temps] = temprep[match(results$dset2[results$dset2 %in% temps], temps)]
      dsets <- param[[which(modal == names(modal_plots))]]$dsets
      results$dset1 <- ordered(results$dset1, levels=names(svals))
      results$dset2 <- ordered(results$dset2, levels=rev(names(svals)))
    } else {
      param$dataset <- as.character(param$dataset)
      param$dataset[param$dataset %in% temps] = temprep[match(param$dataset[param$dataset %in% temps], temps)]
      results <- param[param$modal == modal,]
      results$dataset <- ordered(results$dataset, levels=names(svals))
    }
    if (pp == 1) {
      results$p <- (results$p - minpos[i])/(maxpos[i] - minpos[i])
      results.means <- aggregate(p ~ dataset + modal, data = data.table(results), FUN = mean)
      order <- sort(results.means$p, index.return=TRUE)$ix
      results$dataset <- factor(results$dataset, levels=results.means$dataset[order], ordered=TRUE)
      plot <- ggplot(results, aes(x=p, y=dataset, fill=dataset)) +
        # geom_point(data=results, aes(x=p, y=0, size=size/20, color=dataset), alpha=0.8) +
        geom_density_ridges2(alpha=0.5) +
        xlim(c(minp[i], maxp[i])) +
        scale_fill_manual(name="Study", values=cols) +
        #scale_size_continuous(name="Study", labels=names(svals), breaks=svals) +
        xlab(sprintf("Normalized %s Strength", xl[i])) +
        ylab(TeX("Study")) +
        theme_bw() +
        ggtitle(sprintf("(%s.%s) Study-Wise Distributions", letters[i], roman[1])) +
        labs(shape="Modality", color="Study") +
        theme_bw() +
        theme(legend.position="none", axis.text=element_text(size=11), axis.title=element_text(size=13),
              plot.title = element_text(hjust = -3))
    } else if (files[pp] == 'pair') {
      plot <- ggplot(results, aes(x=dset1, y=dset2, fill=pval)) +
        geom_tile() +
        ggtitle(sprintf(title[pp], letters[i], roman[2])) +
        xlab("Study") +
        ylab("Study") +
        scale_fill_gradientn(name=TeX("$p$-value"), trans="log", breaks=c(0.001, 0.01, 0.1, 1),
                             colours=rev(c("#f2f0f7", "#cbc9e2", "#9e9ac8", "#6a51a3")),
                             limits=c(0.001, 1)) +
        theme_bw() +
        theme(legend.text=element_text(size=10), legend.title=element_text(size=14))
      p_leg <- g_legend(plot)
      plot <- plot +
        theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position="none",
              axis.text=element_text(size=10), axis.title=element_text(size=14),
              plot.title = element_text(hjust = 1.4))
    }
    top_plots[[pp]] <- plot
  }
  top_plot <- grid.arrange(grobs=top_plots, ncol=2)
  cases <- data.frame(case=c(), pval=c(), dataset=c(), size=c())
  for (pp in 3:length(files)) {
    fpath <- paste("data/results", path_str, sprintf('%s.rds', files[pp]), sep="_")
    param = readRDS(fpath)
    results <- param[param$modal == modal,]
    if (pp %in% c(7, 8)) {
      results$dset1 <- as.character(results$dset1)
      results$dset2 <- as.character(results$dset2)
      results$dset1[results$dset1 %in% temps] = temprep[match(results$dset1[results$dset1 %in% temps], temps)]
      results$dset2[results$dset2 %in% temps] = temprep[match(results$dset2[results$dset2 %in% temps], temps)]
      new_dat <- data.frame(case=title[pp], pval=results$pval, dataset=results$dset1, size=results$size)
    } else {
      results$dataset <- as.character(results$dataset)
      results$dataset[results$dataset %in% temps] = temprep[match(results$dataset[results$dataset %in% temps], temps)]
      new_dat <- data.frame(case=title[pp], pval=results$pval, dataset=results$dataset, size=results$size)
    }
    proportion <- rbind(proportion, data.frame(case=title[pp], prop=mean(results$pval < cutoff, na.rm=TRUE)))
    cases <- rbind(cases, new_dat)
  }
  hline = data.frame(y=cutoff, type="sig")
  cases$dataset <- ordered(cases$dataset, levels=names(cols))
  bottom_plot <- ggplot(cases, aes(x=case, y=pval, color=dataset, size=size), drop=FALSE) +
    geom_quasirandom(width=0.25, alpha=0.7) +
    scale_color_manual(name="Study", values=cols) +
    scale_size_continuous(name="Study",  breaks=svals, range=c(0.5, 2)) +
    scale_y_continuous(trans=log10_trans(), limits=c(.001, 1), breaks=c(.001, .01, .1, 1)) +
    #geom_violin(datt=cases, aes(x=case, y=pval, group='black', color='black')) +
    xlab("Pooling Strategy") +
    ylab(TeX("$p$-value")) +
    geom_hline(data=hline, aes(yintercept=y, linetype=type)) +
    scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05"), TeX)) +
    ggtitle(sprintf("(%s.%s) Significance of differences within vs. across each pooling strategy", letters[i], roman[3])) +
    theme_bw() +
    theme(legend.position="none", axis.text=element_text(size=12.5), axis.title=element_text(size=14),
          plot.title = element_text(hjust = -.1), axis.text.x = element_text(size=12.5, angle=20, hjust=1))
  modal_plots[[modal]] <- list(top_plot, bottom_plot)

  case_dat[[modal]] <- cbind(cases, modal)
}

Legends

We make a combined plot just to get the study legend:

case_leg_dat <- do.call(rbind, case_dat)
case_combined <- ggplot(case_leg_dat, aes(x=case, y=pval, color=dataset, size=size)) +
  geom_jitter(alpha=0.8, width=0.25, height=0) +
  scale_color_manual(name="Study", values=cols) +
  scale_size_continuous(name="Study",  breaks=svals, range=c(0.5, 2)) +
  scale_y_continuous(trans=log10_trans(), limits=c(.001, 1), breaks=c(.001, .01, .1, 1)) +
  xlab("Investigation") +
  ylab(TeX("$p$-value")) +
  ggtitle("Batch Effect Investigation") +
  geom_hline(data=hline, aes(yintercept=y, linetype=type)) +
  scale_linetype_manual(values=c("dashed"), name="Cutoff", breaks=c("sig"), labels=lapply(c("$\\alpha = 0.05"), TeX)) +
  theme_bw() +
  theme(legend.text=element_text(size=10), legend.title=element_text(size=14))
case_leg <- g_legend(case_combined)

Combined Plot

dmri <- grid.arrange(grobs=modal_plots$dMRI, heights=c(0.45, 0.55), ncol=1,
                     top=textGrob("dMRI Hemispheric", gp=gpar(fontsize=20)))
fmri <- grid.arrange(grobs=modal_plots$fMRI, heights=c(0.45, 0.55), ncol=1,
                     top=textGrob("fMRI Homotopic", gp=gpar(fontsize=20)))
legs <- grid.arrange(p_leg, case_leg, heights=c(0.45, .55), ncol=1)

grid.arrange(dmri, fmri, legs, widths=c(0.44, 0.44, 0.12))

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