R/post.dist.R

In ANTs: Animal Network Toolkit Software

# Copyright (C) 2018  Sebastian Sosa, Ivan Puga-Gonzalez, Hu Feng He, Xiaohua Xie, Cédric Sueur
#
# This file is part of Animal Network Toolkit Software (ANTs).
#
# ANT is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# ANT is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

#' @title Histogram of posterior distribution
#' @description plot multiple histograms of posterior distribution
#' @param v_perm a data frame of permuted values (each column represents a factor). First row represents the observed values and the next ones represent the permuted values.
#' @param Obs an integer vector indicating the observed value to compare to permuted values if in argument v_perm the first row does not represent the observed values.
#' @param ncols Number of graph per row 
#' @param nrows Number of graph per column
#' @return an histogram.
#' @keywords internal
post.dist <- function(v_perm, Obs = NULL, ncols = NULL, nrows = NULL) {
  par(bg = "gray63")
  if (ncol(v_perm) == 1 & ncol(v_perm) != 0) {
    if (is.null(Obs)) {
      # Extract first values (observed one)
      Obs <- v_perm[1, 1]
      # Extract permuted values
      v_perm <- v_perm[-1, 2]

    }      
    m <- mean(v_perm)
    
    #if (Obs > m) {
    #  # Plot histogram of permuted values
    #  h <- hist(v_perm, breaks = length(v_perm), xaxt = "n")
    #  cuts <- cut(h$breaks, c(Obs, Inf))
    #  cuts <- ifelse(is.na(cuts), "gray10", "gray25")
    #  plot(h, col = cuts, border = cuts, main = paste(attributes(v_perm)$comment), xaxt = "n")
    #  axis(1)
    #  mtext(1, text = round(Obs, digits = 2), at = Obs, col = "white")
    #  # Plot observed value
    #  abline(v = Obs, col = "white")
    #  legend("topright", legend = "observed value", text.col = "white", box.lty = 0)
    #}
    #else {
    #  # Plot histogram of permuted values
    #  h <- hist(v_perm, breaks = length(v_perm), xaxt = "n")
    #  cuts <- cut(h$breaks, c(Obs, Inf))
    #  cuts <- ifelse(is.na(cuts), "gray25", "gray10")
    #  plot(h, col = cuts, border = cuts, xlab = paste(attributes(v_perm)$comment), xaxt = "n")
    #  axis(1)
    #  mtext(1, text = round(Obs, digits = 2), at = Obs, col = "white")
    #  # Plot observed value
    #  abline(v = Obs, col = "white")
    #  legend("topright", legend = "observed value", text.col = "white", box.lty = 0)
    #}
    
    vis.post.distribution(c(Obs, v_perm), legend = FALSE)
  }
  else {
    # Plot output configuration
    if (is.null(ncols) | is.null(nrows)) {
      if (ncol(v_perm) <= 4) {
        (
          par(mfrow = c(1, ncol(v_perm)))
        )
      }
      if (ncol(v_perm) > 4 & ncol(v_perm) < 7) {
        par(mfrow = c(2, 3))
      }
      if (ncol(v_perm) > 6 & ncol(v_perm) < 9) {
        par(mfrow = c(2, 4))
      }
      if(ncol(v_perm) > 9){return(warning("Number of permuted factors are higher than 9. Use ANTs function post.dist to plot them one by one."))}
    }
    else {
      par(mfrow = c(nrows, ncols))
    }
    for (a in 1:ncol(v_perm)) {
      if (is.null(Obs)) {
        # Extract first values (observed one)
        obs <- v_perm[1, a]
        # Extract permuted values
        v_perm <- v_perm[-1, a]
      }
      else{obs=Obs[a]}
      m <- mean(v_perm[, a])
      #if (obs > m) {
      #  # Plot histogram of permuted values
      #  h <- hist(v_perm[, a], breaks = length(v_perm[, a]), plot = FALSE)
      #  cuts <- cut(h$breaks, c(obs, Inf))
      #  cuts <- ifelse(is.na(cuts), "gray10", "gray25")
#
      #  # Plot observed value
      #  if (obs < min(v_perm[, a])) {
      #    xlim <- c(floor(obs / 0.1) * 0.1, ceiling(max(v_perm[, a]) / 0.1) * 0.1)
      #    names(xlim) <- NULL
      #    plot(h, col = cuts, border = cuts, xlim = xlim, main = paste(colnames(v_perm)[a]), xlab = NULL, xaxt = "n")
      #  }
      #  else {
      #    plot(h, col = cuts, border = cuts, main = paste(colnames(v_perm)[a]), xlab = NULL, xaxt = "n")
      #  }
      #  axis(1, pos = 0)
      #  # mtext(1, text = round(obs,digits=3), at = obs, col = "white")
      #  abline(v = obs, col = "white")
      #  legend("topleft", legend = paste("observed value", "\n", round(obs, digits = 3)), text.col = "white", box.lty = 0)
      #}
      #else {
      #  # Plot histogram of permuted values
      #  h <- hist(v_perm[, a], breaks = length(v_perm[, a]), plot = FALSE)
      #  cuts <- cut(h$breaks, c(obs, Inf))
      #  cuts <- ifelse(is.na(cuts), "gray25", "gray10")
      #  # Plot observed value
      #  if (obs > max(v_perm[, a])) {
      #    plot(h, col = cuts, border = cuts, main = paste(colnames(v_perm)[a]), xlab = NULL, xlim = c(floor(v_perm[, a] / 0.1) * 0.1, ceiling(max(obs) / 0.1) * 0.1), xaxt = "n")
      #  }
      #  else {
      #    plot(h, col = cuts, border = cuts, main = paste(colnames(v_perm)[a]), xlab = NULL, xaxt = "n")
      #  }
      #  axis(1, pos = 0)
      #  # mtext(1, text = round(obs,digits=3), at = obs, col = "white")
      #  abline(v = obs, col = "white")
      #  legend("topright", legend = paste("observed value", "\n", round(obs, digits = 3)), text.col = "white", box.lty = 0)
      #}
      vis.post.distribution(c(obs, v_perm[]), legend = FALSE)
    }
  }
}