R/plot.cutter.R

In HelpersMG: Tools for Environmental Analyses, Ecotoxicology and Various R Functions

#' plot.cutter plot result of cutter
#' @title Plot results of cutter that best describe distribution
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Nothing
#' @param x A result file generated by cutter
#' @param col.hist The color of histogram
#' @param col.DL The color of below of above samples
#' @param col.dist The color of distribution
#' @param col.unobserved The color of unobserved states
#' @param col.mcmc The color of mcmc outputs
#' @param legend If TRUE, a legend is shown
#' @param ... Parameters for plot
#' @description Plot the estimates of cut distribution.
#' @family Distributions
#' @examples
#' \dontrun{
#' library(HelpersMG)
#' # _______________________________________________________________
#' # right censored distribution with gamma distribution
#' # _______________________________________________________________
#' # Detection limit
#' DL <- 100
#' # Generate 100 random data from a gamma distribution
#' obc <- rgamma(100, scale=20, shape=2)
#' # remove the data below the detection limit
#' obc[obc>DL] <- +Inf
#' # search for the parameters the best fit these censored data
#' result <- cutter(observations=obc, upper_detection_limit=DL, 
#'                            cut_method="censored")
#' result
#' plot(result, xlim=c(0, 150), breaks=seq(from=0, to=150, by=10))
#' # _______________________________________________________________
#' # The same data seen as truncated data with gamma distribution
#' # _______________________________________________________________
#' obc <- obc[is.finite(obc)]
#' # search for the parameters the best fit these truncated data
#' result <- cutter(observations=obc, upper_detection_limit=DL, 
#'                            cut_method="truncated")
#' result
#' plot(result, xlim=c(0, 150), breaks=seq(from=0, to=150, by=10))
#' # _______________________________________________________________
#' # left censored distribution with gamma distribution
#' # _______________________________________________________________
#' # Detection limit
#' DL <- 10
#' # Generate 100 random data from a gamma distribution
#' obc <- rgamma(100, scale=20, shape=2)
#' # remove the data below the detection limit
#' obc[obc<DL] <- -Inf
#' # search for the parameters the best fit these truncated data
#' result <- cutter(observations=obc, lower_detection_limit=DL, 
#'                           cut_method="censored")
#' result
#' plot(result)
#' plot(result, xlim=c(0, 200), breaks=seq(from=0, to=200, by=10))
#' # _______________________________________________________________
#' # left and right censored distribution
#' # _______________________________________________________________
#' # Generate 100 random data from a gamma distribution
#' obc <- rgamma(100, scale=20, shape=2)
#' # Detection limit
#' LDL <- 10
#' # remove the data below the detection limit
#' obc[obc<LDL] <- -Inf
#' # Detection limit
#' UDL <- 100
#' # remove the data below the detection limit
#' obc[obc>UDL] <- +Inf
#' # search for the parameters the best fit these censored data
#' result <- cutter(observations=obc, lower_detection_limit=LDL, 
#'                            upper_detection_limit=UDL, 
#'                           cut_method="censored")
#' result
#' plot(result, xlim=c(0, 150), col.DL=c("black", "grey"), 
#'                              col.unobserved=c("green", "blue"), 
#'      breaks=seq(from=0, to=150, by=10))
#' # _______________________________________________________________
#' # Example with two values for lower detection limits
#' # corresponding at two different methods of detection for example
#' # with gamma distribution
#' # _______________________________________________________________
#' obc <- rgamma(50, scale=20, shape=2)
#' # Detection limit for sample 1 to 50
#' LDL1 <- 10
#' # remove the data below the detection limit
#' obc[obc<LDL1] <- -Inf
#' obc2 <- rgamma(50, scale=20, shape=2)
#' # Detection limit for sample 1 to 50
#' LDL2 <- 20
#' # remove the data below the detection limit
#' obc2[obc2<LDL2] <- -Inf
#' obc <- c(obc, obc2)
#' # search for the parameters the best fit these censored data
#' result <- cutter(observations=obc, 
#'                            lower_detection_limit=c(rep(LDL1, 50), rep(LDL2, 50)), 
#'                           cut_method="censored")
#' result
#' # It is difficult to choose the best set of colors
#' plot(result, xlim=c(0, 150), col.dist="red", 
#'      col.unobserved=c(rgb(red=1, green=0, blue=0, alpha=0.1), 
#'                       rgb(red=1, green=0, blue=0, alpha=0.2)), 
#'      col.DL=c(rgb(red=0, green=0, blue=1, alpha=0.5), 
#'                       rgb(red=0, green=0, blue=1, alpha=0.9)), 
#'      breaks=seq(from=0, to=200, by=10))
#' # ___________________________________________________________________
#' # left censored distribution comparison of normal, lognormal and gamma
#' # ___________________________________________________________________
#' # Detection limit
#' DL <- 10
#' # Generate 100 random data from a gamma distribution
#' obc <- rgamma(100, scale=20, shape=2)
#' # remove the data below the detection limit
#' obc[obc<DL] <- -Inf
#' # search for the parameters the best fit these truncated data
#' result_gamma <- cutter(observations=obc, lower_detection_limit=DL, 
#'                           cut_method="censored", distribution="gamma")
#' result_gamma
#' plot(result_gamma, xlim=c(0, 250), breaks=seq(from=0, to=250, by=10))
#' 
#' result_lognormal <- cutter(observations=obc, lower_detection_limit=DL, 
#'                           cut_method="censored", distribution="lognormal")
#' result_lognormal
#' plot(result_lognormal, xlim=c(0, 250), breaks=seq(from=0, to=250, by=10))
#' 
#' result_normal <- cutter(observations=obc, lower_detection_limit=DL, 
#'                           cut_method="censored", distribution="normal")
#' result_normal
#' plot(result_normal, xlim=c(0, 250), breaks=seq(from=0, to=250, by=10))
#' 
#' compare_AICc(gamma=result_gamma, 
#'             lognormal=result_lognormal, 
#'             normal=result_normal)
#' # ___________________________________________________________________
#' # Test for similarity in gamma left censored distribution between two
#' # datasets
#' # ___________________________________________________________________
#' obc1 <- rgamma(100, scale=20, shape=2)
#' # Detection limit for sample 1 to 50
#' LDL <- 10
#' # remove the data below the detection limit
#' obc1[obc1<LDL] <- -Inf
#' obc2 <- rgamma(100, scale=10, shape=2)
#' # remove the data below the detection limit
#' obc2[obc2<LDL] <- -Inf
#' # search for the parameters the best fit these censored data
#' result1 <- cutter(observations=obc1, 
#'                   distribution="gamma", 
#'                   lower_detection_limit=LDL, 
#'                   cut_method="censored")
#' logLik(result1)
#' plot(result1, xlim=c(0, 200), 
#'      breaks=seq(from=0, to=200, by=10))
#' result2 <- cutter(observations=obc2, 
#'                   distribution="gamma", 
#'                   lower_detection_limit=LDL, 
#'                   cut_method="censored")
#' logLik(result2)
#' plot(result2, xlim=c(0, 200), 
#'      breaks=seq(from=0, to=200, by=10))
#' result_totl <- cutter(observations=c(obc1, obc2), 
#'                       distribution="gamma", 
#'                       lower_detection_limit=LDL, 
#'                       cut_method="censored")
#' logLik(result_totl)
#' plot(result_totl, xlim=c(0, 200), 
#'      breaks=seq(from=0, to=200, by=10))
#'      
#' compare_AICc(Separate=list(result1, result2), 
#'             Common=result_totl, factor.value=1)
#' compare_BIC(Separate=list(result1, result2), 
#'             Common=result_totl, factor.value=1)           
#' }
#' @method plot cutter
#' @export



plot.cutter <- function(x, 
                        col.hist="grey", 
                        col.DL="blue", 
                        col.dist="black",
                        col.unobserved="green", 
                        col.mcmc=rgb(red=0.6, green = 0.0, blue = 0.0, alpha = 0.01), 
                        legend=TRUE, ...) {
  # x <- result
  # col.hist="grey"
  # col.DL="blue"
  # col.dist="black"
  # col.unobserved="green"
  # legend=TRUE
  # p3p <- list()
  # col.mcmc=rgb(red=0.6, green = 0.0, blue = 0.0, alpha = 0.005)
  
  p3p <- list(...) # p3p <- list()
  
  getparcutter <- getFromNamespace(".getparcutter", ns="HelpersMG")
  
  
  ddistr <- switch(x$distribution, 
                   gamma=dgamma, 
                   lognormal=dlnorm, 
                   normal=dnorm, 
                   weibull=dweibull, 
                   generalized.gamma=dggamma)
  
  qdistr <- switch(x$distribution, 
                   gamma=qgamma, 
                   lognormal=qlnorm, 
                   normal=qnorm, 
                   weibull=qweibull, 
                   generalized.gamma=qggamma)
  
  
  obc <- x$observations[, "Observations"]
  LDL <- unique(x$observations[, "LDL"])
  UDL <- unique(x$observations[, "UDL"])
  
  if (is.null(p3p$xlab))  p3p$xlab <- "Values"
  if (is.null(p3p$ylab))  p3p$ylab <- "Density"
  if (is.null(p3p$main))  p3p$main <- ""
  
  n.mixture <- x$n.mixture
  
  if (!is.null(x$par_mcmc)) {
    par <- x$par_mcmc[2, ]
  } else {
    par <- x$par
  }
  
  parX <- as.list(par)
  pparX <- getparcutter(par, set=NULL)
  parX_mixture <- list()
  for (p in 1:n.mixture)
    parX_mixture <- c(parX_mixture, list(getparcutter(parX, set=p)))
  
  
  if (is.null(p3p$xlim)) {
    nlength <- 100
    mx <- NULL
    for (p in 1:n.mixture) {
      mx <- max(c(mx, do.call(qdistr, args=c(list(p = 0.999, lower.tail = TRUE, log.p = FALSE), parX_mixture[[p]]))))
    }
    
    if (mx > max(x$observations[, "Observations"], na.rm=TRUE) * 5)
      mx <- max(x$observations[, "Observations"], na.rm=TRUE) * 5
    
  } else {
    mx <- p3p$xlim[2]
  }
  
  p3p$xlim <- c(0, mx)
  mi <- 0
  
  if (is.null(p3p$breaks)) {
    es <- suppressWarnings(do.call(hist, args=modifyList(p3p, list(col=col.hist, 
                                                                   x=obc[is.finite(obc)], 
                                                                   plot=FALSE, freq = FALSE))))
    breaks_begin <- es$breaks[1]
    breaks_end <- tail(es$breaks, n=1)
    delta_break <- es$breaks[2]-es$breaks[1]
    
    if (all(!is.na(LDL))) {
      breaks <- unique(c(0, LDL, seq(from=LDL, to=breaks_end+delta_break, by=delta_break)))
    } else {
      breaks <- es$breaks
    }
    if (all(!is.na(UDL))) breaks <- unique(sort(c(UDL, breaks)))
    p3p$breaks <- breaks
  }
  
  npoints <- 10^floor(log10(mx)+1)
  if (npoints < 1000) {
    npoints <- 1000
  }
  x_axis <- seq(from=0, to=mx, length.out=npoints)[-1]
  if (all(!is.na(LDL))) x_axis <- sort(unique(c(LDL, x_axis)))
  if (all(!is.na(UDL))) x_axis <- sort(unique(c(UDL, x_axis)))
  npoints <- length(x_axis) + 1
  
  y_axis <- rep(0, npoints-1)
  for (p in 1:n.mixture)
    y_axis <- y_axis + do.call(ddistr, args=c(list(x=x_axis, log = FALSE), parX_mixture[[p]])) * pparX[p]
  
  if ((!is.null(x$mcmc)) & (!is.null(col.mcmc))) {
    pl_mcmc <- matrix(data = NA, ncol=length(x_axis), nrow = nrow(x$mcmc$resultMCMC[[1]]))
    
    # Là je fais tous les MCMC
    for (i in 1:nrow(x$mcmc$resultMCMC[[1]])) {
      par_mcmc <- x$mcmc$resultMCMC[[1]][i, ]
      pparX_mcmc <- getparcutter(par_mcmc, set=NULL)
      
      y <- rep(0, length(x_axis))
      for (p in 1:n.mixture)
        y <- y + do.call(ddistr, args=c(list(x=x_axis, log = FALSE), as.list(getparcutter(par_mcmc, set=p)))) * pparX_mcmc[p]
      
      pl_mcmc[i, ] <- y
    }
  } else {
    pl_mcmc <- NULL
  }
  
  
  if (is.null(p3p$ylim)) {
    es <- suppressWarnings(do.call(hist, args=modifyList(p3p, list(col=col.hist, 
                                                                   x=obc[is.finite(obc)], plot=FALSE, freq = FALSE))))
    y <- y_axis
    if (!is.null(pl_mcmc)) y <- c(y, as.vector(pl_mcmc))
    
    p3p$ylim <- c(0, max(c(y[is.finite(y)], es$density), na.rm=TRUE))
  }
  
  es <- do.call(hist, args=modifyList(p3p, list(col=col.hist, 
                                                x=obc[is.finite(obc)], plot=TRUE, freq = FALSE)))
  
  
  # Dans u j'ai le nombre total de limites
  u <- 0
  if (all(!is.na(UDL))) u <- length(UDL)
  if (all(!is.na(LDL))) u <- u + length(LDL)
  
  if (u == 0) {
    if (legend) {
      if (is.null(pl_mcmc)) {
        legend("topright", legend = c("Observed values", "Fitted distribution"), 
               col=c(col.hist, col.dist), 
               pch=c(15, NA), lty=c(NA, 1))
      } else {
        c2 <- col.mcmc
        substr(c2, 8, 9) <- "88"
        legend("topright", legend = c("Observed values", "Fitted distribution", "Posterior predictive distribution"), 
               col=c(col.hist, col.dist, c2), 
               pch=c(15, NA, NA), lty=c(NA, 1, 1), lwd=c(NA, 1, 10))
      }
    }
  } else {
    
    col.unobserved_X <- rep(col.unobserved, u)[1:u]
    col.DL_X <- rep(col.DL, u)[1:u]
    cpt_u <- 0
    
    if (all(!is.na(LDL))) {
      
      for (i in seq_along(LDL)) {
        posy <- i/(length(LDL)+1)
        x_axis_LDL <- x_axis[x_axis<=LDL[i]]
        # npoints <- 100
        # x_axis_LDL <- seq(from=0, to=LDL[i], length.out=npoints)
        y_LDL <- rep(0, length(x_axis_LDL))
        
        for (p in 1:n.mixture)
          y_LDL <- y_LDL + do.call(ddistr, args = modifyList(list(x=x_axis_LDL, log = FALSE), 
                                                             parX_mixture[[p]])) * pparX[p]
        
        # y_LDL <- ifelse(y_LDL>max(y), max(y), y_LDL)
        polygon(c(x_axis_LDL, LDL[i], 0), 
                c(y_LDL, 0, 0), 
                col=col.unobserved_X[cpt_u+i], border=col.unobserved_X[cpt_u+i])
        
      }
      
      if (!is.null(x$LDL)) {
        for (i in seq_along(LDL)) {
          posy <- i/(length(LDL)+1)
          
          if (nrow(x$LDL) == 3) {
          
          segments(x0=x$LDL[1, i], x1=x$LDL[3, i], y0=ScalePreviousPlot(x=0, y=posy)$y, 
                   y1=ScalePreviousPlot(x=0, y=posy)$y, 
                   col=col.DL_X[cpt_u+i])
          points(x=x$LDL[2, i], y=ScalePreviousPlot(x=0, y=posy)$y, pch=19, col=col.DL_X[cpt_u+i])
          segments(x0=x$LDL[1, i], x1=x$LDL[1, i], y0=ScalePreviousPlot(x=0, y=posy-0.03)$y, y1=ScalePreviousPlot(x=0, y=posy+0.03)$y, 
                   col=col.DL_X[cpt_u+i])
          segments(x0=x$LDL[3, i], x1=x$LDL[3, i], y0=ScalePreviousPlot(x=0, y=posy-0.03)$y, y1=ScalePreviousPlot(x=0, y=posy+0.03)$y, 
                   col=col.DL_X[cpt_u+i])
          } else {
            points(x=x$LDL[1, i], y=ScalePreviousPlot(x=0, y=posy)$y, pch=19, col=col.DL_X[cpt_u+i])
          }
        }
        cpt_u <- cpt_u + length(LDL)
      }
    }
    
    if (all(!is.na(UDL))) {
      for (i in seq_along(UDL)) {
        posy <- i/(length(UDL)+1)
        npoints <- 100
        x_axis_UDL <- seq(from=UDL[i], to=mx, length.out=npoints)
        
        y_UDL <- rep(0, length(x_axis_UDL))
        for (p in 1:n.mixture)
          y_UDL <- y_UDL + do.call(ddistr, args = modifyList(list(x=x_axis_UDL, log = FALSE), 
                                                             parX_mixture[[p]])) * pparX[p]
        
        # y_UDL <- ifelse(y_UDL > max(y), max(y), y_UDL)
        polygon(c(x_axis_UDL, mx, UDL[i]), 
                c(y_UDL, 0, 0), 
                col=col.unobserved_X[cpt_u+i], border=col.unobserved_X[cpt_u+i])
      }
      
      if (!is.null(x$UDL)) {
        for (i in seq_along(UDL)) {
          posy <- i/(length(UDL)+1)
          
          if (nrow(x$UDL) == 3) {
          segments(x0=x$UDL[1, i], x1=x$UDL[3, i], y0=ScalePreviousPlot(x=0, y=posy)$y, y1=ScalePreviousPlot(x=0, y=posy)$y, 
                   col=col.DL_X[cpt_u+i])
          points(x=x$UDL[2, i], y=ScalePreviousPlot(x=0, y=posy)$y, pch=19, col=col.DL_X[cpt_u+i])
          segments(x0=x$UDL[1, i], x1=x$UDL[1, i], y0=ScalePreviousPlot(x=0, y=posy-0.03)$y, y1=ScalePreviousPlot(x=0, y=posy+0.03)$y, 
                   col=col.DL_X[cpt_u+i])
          segments(x0=x$UDL[3, i], x1=x$UDL[3, i], y0=ScalePreviousPlot(x=0, y=posy-0.03)$y, y1=ScalePreviousPlot(x=0, y=posy+0.03)$y, 
                   col=col.DL_X[cpt_u+i])
          } else {
            points(x=x$UDL[1, i], y=ScalePreviousPlot(x=0, y=posy)$y, pch=19, col=col.DL_X[cpt_u+i])
          }
        }
      }
    }
    
    
    if (legend) {
      if (is.null(pl_mcmc)) {
        legend("topright", legend = c("Observed values", "Fitted distribution", 
                                      paste(rep("Unobserved distribution", u), na.omit(c(LDL, UDL))), 
                                      paste(rep("Mean unobserved", u), na.omit(c(LDL, UDL)))), 
               col=c(col.hist, col.dist, col.unobserved_X, col.DL_X), 
               pch=c(15, NA, rep(17, u), rep(19, u)), lty=c(NA, 1, rep(NA, u) , rep(1, u)))
      } else {
        c2 <- col.mcmc
        substr(c2, 8, 9) <- "88"
        legend("topright", legend = c("Observed values", "Fitted distribution", "Posterior predictive distribution", 
                                      paste(rep("Unobserved distribution", u), na.omit(c(LDL, UDL))), 
                                      paste(rep("Mean unobserved", u), na.omit(c(LDL, UDL)))), 
               col=c(col.hist, col.dist, c2, col.unobserved_X, col.DL_X), 
               pch=c(15, NA, NA, rep(17, u), rep(19, u)), lty=c(NA, 1, 1, rep(NA, u) , rep(1, u)), 
               lwd=c(NA, 1, 10, rep(NA, u) , rep(1, u)))
        
      }
    }
  }
  
  if (!is.null(pl_mcmc)) {
    for (i in 1:nrow(pl_mcmc)) {
      y <- pl_mcmc[i, ]
      lines(x_axis[is.finite(y)], y[is.finite(y)], col=col.mcmc)
    }
  }
  
  lines(x_axis[is.finite(y_axis)], y_axis[is.finite(y_axis)], col=col.dist)
  
}