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)
source('../../R/siem.R')

localpath <- '/home/eric/Documents/research/R-fmri/graphstats'


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('colors.rds')
svals <- colors$svals; cols <- colors$cols

Aggregate plots

panel_plots <- list()
files <- c("params", "pair", "case1", "case2", "case3", "case4", "case5", "case6")
title <- c("Exploratory", "Paired", "1", "2", "3", "4",
           "5", "6")
modal_plots <- list(dMRI=list(), fMRI=list())
counters <- list(dMRI=1, fMRI=1)
for (pp in 1:length(files)) {
  cases <- data.frame(case=c(), dataset=c(), pval=c(), size=c())
  for (modal in names(modal_plots)) {
    if (modal == 'dMRI') {
      path_str <- 'hem'
    } else {
      path_str <- 'bil'
    }
    fpath <- file.path(localpath, 'docs/siem',
                       paste("results", path_str, sprintf('%s.rds', files[pp]), sep="_"))
    param = readRDS(fpath)
    if (pp == 2) {
      results <- param[[which(modal == names(modal_plots))]]$data
      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 {
      results <- param[param$modal == modal,]
      if (pp %in% c(7, 8)) {
        results$dset1 <- ordered(results$dset1, levels=names(svals))
        results$dset2 <- ordered(results$dset2, levels=names(svals))
      } else {
        results$dataset <- ordered(results$dataset, levels=names(svals))
      }
    }
    if (pp == 2) {
      plot <- ggplot(results, aes(x=p, y=q, size=size, color=dataset)) +
        geom_point(alpha=0.8) +
        scale_color_manual(name="Dataset", values=cols) +
        scale_size_continuous(name="Dataset", labels=names(svals), breaks=svals) +
        xlab(TeX("$\\hat{p}$")) +
        ylab(TeX("$\\hat{q}$")) +
        theme_bw() +
        ggtitle("Exploratory Analysis") +
        labs(shape="Modality", color="Dataset")
    } else if (files[pp] == 'pair') {
      plot <- ggplot(results, aes(x=dset1, y=dset2, fill=pval)) +
        geom_tile() +
        ggtitle(title[pp]) +
        xlab("Dataset") +
        ylab("Dataset") +
        scale_fill_gradientn(name=TeX("$p$-value"), trans="log", breaks=c(0.001, 0.01, 0.1, 1),
                             colours=c("#f2f0f7", "#cbc9e2", "#9e9ac8", "#6a51a3"), limits=c(0.001, 1)) +
        theme(axis.text.x = element_text(angle = 90, hjust = 1))
    } else if (files[pp] == 'case2') {
      plot <- ggplot(results, aes(pval, color=dataset)) +
        geom_density(alpha=0.2, size=1.5) +
        scale_color_manual(name="Dataset", values=cols) +
        xlab("") +
        coord_flip() +
        ylab("") +
        theme_bw() +
        scale_x_continuous(trans=log10_trans(), limits=c(.001, 1), breaks=c(.001, .001, .01, .1, 1)) +
        ggtitle(title[pp])
    }  else {
      if (pp %in% c(7, 8)) {
        plot <- ggplot(results, aes(x=modal, y=pval, color=dset1, size=size)) +
          scale_size_continuous(name="Dataset", breaks=svals)
      } else if (!(pp %in% c(6))) {
        plot <- ggplot(results, aes(x=modal, y=pval, color=dataset, size=size)) +
          scale_size_continuous(name="Dataset", breaks=svals)
      } else {
        plot <- ggplot(results, aes(x=modal, y=pval, color=dataset))
      }
      plot <- plot +
        geom_jitter(alpha=0.8, width=0.25, height=0) +
        scale_color_manual(name="Dataset", values=cols) +
        ggtitle(title[pp]) +
        ylab(TeX("$p$-Value")) +
        xlab("Modality") +
        scale_y_continuous(trans=log10_trans(), limits=c(.001, 1), breaks=c(.001, .001, .01, .1, 1)) +
        labs(shape="Modality", color="Dataset") +
        theme_bw() +
        ggtitle(title[pp])
    }
    modal_plots[[modal]][[pp]] <- plot
  }
}

We make a combined plot just to get the legend:

leg_dat <- readRDS('./results_bil_params.rds')
leg_dat$dataset <- ordered(leg_dat$dataset, levels=names(svals))
leg_plot <- ggplot(leg_dat, aes(x=p, y=q, size=size, color=dataset)) +
  geom_point(alpha=1) +
  scale_color_manual(name="Dataset", values=cols) +
  scale_size_continuous(name="Dataset", breaks=svals) +
  xlab(TeX("$\\hat{p}$")) +
  ylab(TeX("$\\hat{q}$")) +
  theme_bw() +
  ggtitle("Exploratory Analysis") +
  labs(shape="Modality", color="Dataset") +
  guides(col=guide_legend(ncol=1))
leg_ex <- g_legend(leg_plot)

Combined Plot

comb_plot <- list()
title <- c("dMRI, Hemispheric", "fMRI, Homotopic")
for (i in 1:length(modal_plots)) {
  modal <- modal_plots[[i]]
  modal <- lapply(1:length(modal), function(j) {
      plot <- modal[[j]]
      plot <- plot + theme(legend.position="none")
      if ((j >= 5) & !(j %in% c(4))) {
        plot <- plot + xlab("") + ylab("")
      }
      if (j > 2) {
        plot <- plot + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
      }
      return(plot)
    })
  comb_plot[[i]] <- grid.arrange(grid.arrange(grid.arrange(grobs=list(modal[[1]], modal[[2]]), nrow=1),
                            grid.arrange(grobs=modal[c(3, 4, 5, 6, 7, 8)], nrow=1), heights=c(0.5, 0.5)),
                            top=textGrob(title[i], gp=gpar(fontsize=20)))
}

grid.arrange(grobs=c(comb_plot, list(leg_ex)), ncol=3, widths=c(0.44, 0.44, 0.12))

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