R/fmri_pval_comparison_2d.R

In TCIU: Spacekime Analytics, Time Complexity and Inferential Uncertainty

#' @title 2D comparison visualization between the p-values
#' @description a plot arrangement method, which uses \code{gridExtra} to combine multiple 
#' 2D plots of the fMRI data together. This can bring convenience for users to compare 
#' the result of different statistical tests based on the p values they provide
#'
#' @param pval_ls a list. Each element is a 3D array of p-values data
#' @param pval_name_ls a list with the element as name for the p-values data provided in 'pval_ls'
#' @param axis_i_lses a list with 3 numeric elements or a list of lists. If the elements are numeric, 
#' they would specify indices of slice for the three direction. 
#' If any direction of the slice need not to be shown, make it as NULL for that element. 
#' If elements are lists, each list provides specified cuts for corresponding 3D p-values data.
#' @param hemody_data a parameter to have the plot with/without hemodynamic contour. The default is NULL to make the plot 
#' without hemodynamic contour, otherwise assign a 3D array of the hemodynamic data.
#' @param mask a 3D nifti or 3D array of data to show the shell of the brain
#' @param p_threshold NULL or a numeric value that can be selected randomly below 0.05 to 
#' drop all p-values above the threshold. If 'low5_percent' method is used, 
#' make 'p_threshold' as NULL. The default is 0.05.
#' @param legend_show a logical parameter to specify whether the final plot has legend
#' for all the subplots or the shared legend for all the subplots. The default is TRUE.
#' @param method a string that represents method for the plot. 
#' There are 3 options: 'min_max', 'scale_p' and 'low5_percent'. The default is 'scale_p'. 
#' 'min_max' is to draw plot based on the color scale of the minimum and maximum of the p-values; 
#' 'scale_p' is to draw the plot with fixed color scale for fixed range of p-values; 
#' 'low5_percent' is to draw the plot for the smallest 5 percent of p-values
#' when all the p-values are not significant
#' @param color_pal the name of the color palettes provided by \code{RColorBrewer} The default is "YlOrRd".
#' @param multi_pranges an option under 'scale_p' method to decide whether there are at most 9 colors 
#' in the legend for the ranges of p-values, or at most 4 colors. 
#' The default is TRUE, choosing the larger number of colors for the plot.
#' @param mask_width a numeric value to specify the width of mask contour. The default is 1.5.
#' 
#' @details
#' The function \code{fmri_pval_comparison_2d} is used to combine and compare the 2D plots for different 3D arrays of p-values.
#' The plots in each row are generated by one specific 3D p value data. 
#' The first column of the integrated plot specifies the name of the 3D p value data (for generation of 
#' the plots in that row). The rest of the three columns are the plots from sagittal, coronal and 
#' axial view for each 3D p value data.
#'
#' @author SOCR team <\url{http://socr.umich.edu/people/}>
#'
#' @return a combination plot arranged by \code{gridExtra}
#'
#' @examples
#' # sample 3D data of mask provided by the package
#' dim(mask)
#' # sample 3D p value provided by the package
#' dim(phase2_pval)
#' dim(phase3_pval)
#' 
#' \donttest{
#' fmri_pval_comparison_2d(list(phase2_pval, phase3_pval), 
#'                         list('phase2_pval', 'phase3_pval'),
#'                         list(list(40, 26, 33), list(40, 26, 33)), 
#'                         hemody_data = NULL, 
#'                         mask = mask, p_threshold = 0.05, 
#'                         legend_show = FALSE, method = 'scale_p',
#'                         color_pal = "YlOrRd", multi_pranges=TRUE)
#' }                  
#' @export
#'
#' @import dplyr scales RColorBrewer ggpubr
#' @importFrom ggplot2 ggplot geom_line scale_fill_identity labs geom_contour ggtitle theme geom_tile scale_alpha scale_fill_gradient coord_fixed aes
#' @importFrom gridExtra grid.arrange arrangeGrob

fmri_pval_comparison_2d = function(pval_ls,
                                   pval_name_ls,
                                   axis_i_lses,
                                   hemody_data = NULL, 
                                   mask,
                                   p_threshold = 0.05, 
                                   legend_show = TRUE,
                                   method = "scale_p",
                                   color_pal = "YlOrRd",
                                   multi_pranges = TRUE,
                                   mask_width = 1.5){
  
  if((inherits(axis_i_lses,"list")) &
     (inherits(axis_i_lses[[1]],"list"))){
    axis_i_ls = do.call(c, axis_i_lses)
    if(length(axis_i_ls)!=3*length(axis_i_lses)){
      stop("If 'axis_i_ls'contains more than one list of index",
           ", ", "each list should have 3 elements for x, y, z.")
    }
  }else if((inherits(axis_i_lses,"list")) &
           (inherits(axis_i_lses[[1]],"numeric")) &#change, because the first element can be null
           (length(axis_i_lses) == 3)){
    axis_i_ls = rep(axis_i_lses, length(pval_ls))
  }else if((inherits(axis_i_lses,"list") == FALSE)| 
           (length(axis_i_lses)!=3)){
    stop("'axis_i_ls' should be list with 3 elements for x, y, z.")
  }###!!! not accurate, change!
  
  # generate cutomized error message based on different wrong input
  if((inherits(pval_ls,"list")==FALSE) |
     (length(pval_ls)!=length(pval_name_ls)) |
     (inherits(pval_name_ls,"list")==FALSE)
  ){
    stop("'pval_ls', 'pval_name_ls'should be list with the same number of  elements.")
  }else if(is.null(hemody_data) != TRUE){
    if((inherits(hemody_data,"array")==FALSE) | 
       (length(dim(hemody_data))!=3)){
      stop("If 'hemody_data' is not NULL, then it should be a 3d array.")
    }
  }else if(((class(mask) %in% c("array", "nifti")) != TRUE) |
           (length(dim(mask))!=3)){
    stop("'mask' should be a 3d array.")
  }else if((inherits(p_threshold,"numeric") ==FALSE) | 
           (p_threshold > 0.05) | (p_threshold <= 0)){
    stop("'p_threshold should be a numeric value in range of (0, 0.5].'")
  }else if(inherits(legend_show,"logical") ==FALSE){
    stop("'legend_show' should be a logical TRUE or FALSE.")
  }else if((method %in% c("scale_p", "min_max")) != TRUE){
    stop("'method' should only choose from 'scale_p' or 'min_max'.")
  }
  
  len_idx_sets = length(axis_i_ls)/3
  # get the specific numeric locations for x, y, z respectively to slice the 3d fMRI data  
  axis_vec = unique(rep(c("x", "y", "z"), 
                    len_idx_sets)[unlist( lapply( axis_i_ls, 
					                              function(i) ifelse(is.null(i),FALSE,TRUE)) )])
  axis_i_vec = unlist(axis_i_ls[unlist( lapply( axis_i_ls, 
                                                function(i) ifelse(is.null(i),FALSE,TRUE)) )])
  
  # get the plots for every 3d p values provided from every view that specified
  pl_name_vec = c()
  ing_pls_name_vec = c()
  idx = 0
  for( i in c(1:length(pval_ls)) ){
    
    for( axis_ele in axis_vec ){
      idx = (idx + 1)
      # print(idx)
      # give error when the 2d plot function does not work successfully
      pl = tryCatch({
        fmri_2dvisual(pval_ls[[i]], list(axis_ele, axis_i_vec[idx]), 
                      hemody_data, mask, p_threshold, legend_show,
                      method, color_pal, multi_pranges, mask_width)},
        error = function(cond){
          message("The mistake is at the 2d plot function.")
          message(paste("The message error in 2d plot is: ", cond))
          message(paste("The wrong plot is "), paste0(axis_ele, as.character(i)))
          message("The number represents the plot for which p values", 
		          ", ", 
				  "and x, y, z represents sagittal, coronal and axial view.")
          # return(NA)
          return( ggpubr::text_grob(as.character("NULL")) )
        }
      )
      pl_name = paste0(axis_ele, as.character(i))
      # assign the 2d plot to the specified name whose type is character
      assign(pl_name, pl) 
      
      pl_name_vec = append(pl_name_vec, pl_name)
    }
    
    ing_pls_name = paste0("pls", as.character(i))
    ing_pls_name_vec = append(ing_pls_name_vec, ing_pls_name)
  }
  
  
  # combine the plots from the same 3d p value data in one row in one plot
  idx=0
  
  for (plt_item in ing_pls_name_vec){
    start_idx = idx*length(axis_vec)+1
    end_idx = (idx+1)*length(axis_vec)
    int_plt_str = pl_name_vec[start_idx: end_idx]
    idx = idx + 1
    text_str = pval_name_ls[[idx]]
    
    int_plt_str_trueval = sapply(int_plt_str, get, envir=sys.frame(sys.parent(0)), simplify=FALSE)
    int_plts_form = append( int_plt_str_trueval,
                            list(nrow=1, ncol=1+length(axis_vec)) )
    int_plts_form_wtext = append( list(ggpubr::text_grob(text_str)), int_plts_form )
    int_plt = do.call("arrangeGrob", int_plts_form_wtext)
    assign(plt_item, int_plt)
  }
  
  
  
  if(legend_show == FALSE){
    
    # below is to get the shared legend
    # break points to generate intervals of 1 - p value
    complement_p_cut = c(1-5e-2, 1-1e-2, 1-1e-3, 1-1e-4, 1-1e-5,#1-1e-0, 
                         1-1e-6, 1-1e-7, 1-1e-8, 1-0)
    # get all possible 1-p value interval
    one_minus_p_range = levels(cut(c(0, 1), breaks = complement_p_cut, right = FALSE))
    # get p intervals base on 1-p intervals
    numeric_range_one_minus_p = regmatches(one_minus_p_range, 
                                           gregexpr("[[:digit:]\\.]+", 
                                                    one_minus_p_range))
    label_p_range = c()
    for(i in 1:length(numeric_range_one_minus_p)){
      # formatC(1-numeric_range_one_minus_p[[i]], format = "e", digits = 1)
      p_2range = 1 - as.numeric(numeric_range_one_minus_p[[i]])
      label_p_range[i] = sprintf("(%.1e, %.1e]", p_2range[2], p_2range[1])
    }
    # decide the color for each voxel based on its 1-p range interval
    p_range_num = as.numeric(sort(as.factor(label_p_range)))
    one_minus_p_color = brewer.pal(9,"YlOrRd")[p_range_num]
    p_range_num_df = data.frame(x = p_range_num, y = p_range_num, 
                                color_i = one_minus_p_color)
    # make the plot to extract its legend
    color_i = one_minus_p_color
    p = ggplot() +
      geom_tile(p_range_num_df, mapping=aes(x=x, y=y, fill=color_i)) + 
      scale_fill_identity("p value", labels= label_p_range, 
                          breaks = sort(unique(p_range_num_df$color_i), 
                                        decreasing = TRUE), 
                          guide = "legend") + theme(legend.position = "bottom")
    legend_p = get_legend(p) 
    
    # integrate all plots (every row of plots) together with one shared legend
    ing_pls_name_trueval = sapply(ing_pls_name_vec, get, 
	                              envir=sys.frame(sys.parent(0)), simplify=FALSE)
    ing_pls_name_form = append(ing_pls_name_trueval, list(legend_p))
    plts_total = do.call("grid.arrange", append( ing_pls_name_form,
                                                 list(nrow=(1+length(pval_ls))) )) + coord_fixed()
  }else{
    
    # combine all the subplots with their own legend together
    ing_pls_name_trueval = sapply(ing_pls_name_vec, get, 
	                              envir=sys.frame(sys.parent(0)), simplify=FALSE)
    plts_total = do.call("grid.arrange", append( ing_pls_name_trueval,
                                                 list(nrow=length(pval_ls)) )) + coord_fixed()
  }
  
}