R/print.cutter.R

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

#' print.cutter plot result of cutter
#' @title Print 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 silent If TRUE does not show the outpout
#' @param ... Not used
#' @description Print 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, 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 print cutter
#' @export



print.cutter <- function(x, silent=FALSE, ...) {
  
  # p3p <- list(...)
  
  mc <- "censored and truncated"
  if (all(x$observation[, "Cut"] == "censored")) mc <- "censored"
  if (all(x$observation[, "Cut"] == "truncated")) mc <- "truncated"
  
  textoutput <- paste(x$distribution, mc, "distribution using Bayesian MCMC:")
  textoutput <- paste0(textoutput, "\n", paste0(rep("-", nchar(textoutput)), collapse = ""))
  
  textoutput <- paste0(textoutput, "\nNumber of data: ", x$n.quantified+x$n.LDL+x$n.UDL)
  textoutput <- paste0(textoutput, "\nNumber of LDL data: ", x$n.LDL, " (", specify_decimal(x$n.LDL/(x$n.quantified+x$n.LDL+x$n.UDL)), ")")
  textoutput <- paste0(textoutput, "\nNumber of UDL data: ", x$n.UDL, " (", specify_decimal(x$n.UDL/(x$n.quantified+x$n.LDL+x$n.UDL)), ")")
  if (!is.null(x$prob.LDL)) {
    textoutput <- paste0(textoutput, "\nExpected proportion of LDL data: ")
    textoutput <- paste0(textoutput, "\n", paste0(capture.output(print(x$prob.LDL)), collapse = "\n"))
  }
  if (!is.null(x$prob.UDL))  {
    textoutput <- paste0(textoutput, "\nExpected proportion of UDL data: ")
    textoutput <- paste0(textoutput, "\n", paste0(capture.output(print(x$prob.UDL)), collapse = "\n"))
  }
  
  if (!is.null(x$par_mcmc)) {
    for (i in 1:ncol(x$par_mcmc)) {
      textoutput <- paste0(textoutput, "\n", paste(colnames(x$par_mcmc)[i], "=", specify_decimal(x$par_mcmc[2, i]), "; Credible Interval 95%= [", specify_decimal(x$par_mcmc[1, i]), ",", specify_decimal(unname(x$par_mcmc[3, i])), "]", sep=""))
    }
  } else {
    for (i in 1:length(x$par)) {
      textoutput <- paste0(textoutput, "\n", paste(names(x$par)[i], "=", specify_decimal(x$par[i]), sep=""))
    }
  }
  
  
  textoutput <- paste0(textoutput, "\n", paste("-Ln L=", specify_decimal(x$value)))
  textoutput <- paste0(textoutput, "\n", paste("AIC=", specify_decimal(x$AIC)))
  textoutput <- paste0(textoutput, "\n", paste("AICc=", specify_decimal(x$AICc)))
  textoutput <- paste0(textoutput, "\n", paste("BIC=", specify_decimal(x$BIC)))
  
  mean <- x$total.distribution["mean"]
  sd<- x$total.distribution["sd"]
  
  textoutput <- paste0(textoutput, "\n", paste0("The distribution has a mean of ", specify_decimal(mean), 
                                                " and a standard deviation of ", specify_decimal(sd), "."))
  
  if (!is.null(x$LDL)) {
    for (LDL_c in 1:ncol(x$LDL)) {
      if (nrow(x$LDL)==3) {
      textoutput <- paste0(textoutput, "\n", paste0("The median value for samples below the detection limit ", specify_decimal(colnames(x$LDL)[LDL_c])," is ", specify_decimal(x$LDL[2, LDL_c])))
      textoutput <- paste0(textoutput, "\n", paste0("   The 95% credible interval is ", specify_decimal(x$LDL[1, LDL_c]), " to ", specify_decimal(x$LDL[3, LDL_c])))
      } else {
        textoutput <- paste0(textoutput, "\n", paste0("The mean value for samples below the detection limit ", specify_decimal(colnames(x$LDL)[LDL_c])," is ", specify_decimal(x$LDL["mean", LDL_c])))
      }
    }
  }
  
  if (!is.null(x$UDL)) {
    for (UDL_c in 1:ncol(x$UDL)) {
      if (nrow(x$UDL)==3) {
      textoutput <- paste0(textoutput, "\n", paste0("The median value for samples above the detection limit ", specify_decimal(colnames(x$UDL)[UDL_c])," is ", specify_decimal(x$UDL[2, UDL_c])))
      textoutput <- paste0(textoutput, "\n", paste0("   The 95% credible interval is ", specify_decimal(x$UDL[1, UDL_c]), " to ", specify_decimal(x$UDL[3, UDL_c])))
      } else {
        textoutput <- paste0(textoutput, "\n", paste0("The mean value for samples above the detection limit ", specify_decimal(colnames(x$UDL)[UDL_c])," is ", specify_decimal(x$UDL["mean", UDL_c])))
      }
    }
  }
  
  if (!is.null(x$GoF)) {
    
    textoutput <- paste0(textoutput, "\n", paste("The probability that the data could have been obtained using the fitted distribution is", specify_decimal(x$GoF), "taking into account the distribution of likelihoods of posterior predictive distribution compared to the likelihood of observed data in the fitted model."))
    
  }
  textoutput <- paste0(textoutput, "\n")
  
  if (!silent) cat(textoutput)
  return(invisible(textoutput))
}