R/plot.ECFOCF.R

In phenology: Tools to Manage a Parametric Function that Describes Phenology and More

#' plot.ECFOCF plots a result of clutch frequency fit.
#' @title Plot a result of clutch frequency fit.
#' @author Marc Girondot
#' @return Nothing
#' @param x A result for fitCF().
#' @param ... Graphic parameters, see plot.TableECFOCF() or par.
#' @param result What result will be plotted: data, dataOCF, dataECF, ECF, OCF, ECFOCF, ECFOCF0, CF, Probabilities, period
#' @param parameters Name of parameters to plot for Probabilities result.
#' @param legend Legend to use for the parameters as a vector of names.
#' @param category What category will be plotted, numeric or NA for all.
#' @param period The period that will be plotted.
#' @param resultMCMC A result from fitRMU_MHmcmc.
#' @param chain Which chain to be used in resultMCMC.
#' @param replicates How many replicates fron resultMCMC.
#' @description This function plots the result of fitCF().\cr
#' The result \code{data} plots the observed ECF-OCF table.\cr
#' The result \code{dataOCF} plots the observed OCF table.\cr
#' The result \code{dataECF} plots the observed ECF table.\cr
#' The result \code{CF} plots the true clutch frequency.\cr
#' The result \code{OCF} plots the observed clutch frequency.\cr
#' The result \code{ECF} plots the estimated clutch frequency.\cr
#' The result \code{ECFOCF} plots the bivariate observed vs. estimated clutch frequency.\cr
#' The result \code{ECFOCF0} plots the bivariate observed vs. estimated clutch frequency without the 0 OCF.\cr
#' The result \code{probabilities} plots the probabilities of captures, p and a, or OTN.\cr
#' The result \code{period} plots the probabilities of nesting according to period.\cr
#' If category is left to NA, the compound value for all the population is plotted.\cr
#' When result="data" is used, this is a parser for plot.TableECFOCF().\cr
#' See this function for the parameters.\cr
#' The parameter y.axis is the shift of the x legends for result="prob".\cr
#' When a \code{resultMCMC} is indicated, if replicates is "all", all values from all chains are used; 
#' if a value lower than number of iterations is indicated, a regular thinning is used and 
#' if a value larger then number if iteration is indicated, a sampling with replacement is used.\cr
#' @family Model of Clutch Frequency
#' @examples
#' \dontrun{
#' library(phenology)
#' # Example
#' data(MarineTurtles_2002)
#' ECFOCF_2002 <- TableECFOCF(MarineTurtles_2002)
#' o_mu1p2_NB <- fitCF(x = c(mu = 4.6426989650675701, 
#'                          sd = 75.828239144717074, 
#'                          p1 = 0.62036295627161053,
#'                          p2 = -2.3923021862881511, 
#'                          OTN = 0.33107456308054345),
#'                  data=ECFOCF_2002)
#'                  
#' par(mar=c(4, 4, 1, 1)+0.4)
#' plot(o_mu1p2_NB, result="data", category=NA, 
#'      bty="n", las=1, cex.points=3, cex.axis = 0.8)
#' plot(o_mu1p2_NB,result="data", category=NA, 
#'      bty="n", las=1, cex.points=3, pch=NA, 
#'      col.labels = "red", show.labels=TRUE, cex.0=0.2, 
#'      show.0 = TRUE, col.0="blue", pch.0=4)
#' plot(o_mu1p2_NB, result="dataOCF", category=NA, 
#'      bty="n", las=1)
#' plot(o_mu1p2_NB, result="dataECF", category=NA, 
#'      bty="n", las=1)
#'      
#' plot(o_mu1p2_NB, result="CF", bty="n", las=1)
#' 
#' plot(o_mu1p2_NB, result="OCF", category=1, bty="n", las=1)
#' plot(o_mu1p2_NB, result="OCF", category=2, bty="n", las=1)
#' 
#' plot(o_mu1p2_NB, result="ECFOCF", bty="n", las=1)
#' 
#' plot(o_mu1p2_NB, result="ECFOCF0", bty="n", las=1)
#' plot(o_mu1p2_NB, result="ECFOCF0", category=1, bty="n", las=1)
#' plot(o_mu1p2_NB, result="ECFOCF0", category=2, bty="n", las=1)
#' 
#' plot(o_mu1p2_NB, result="Prob", category=c(1, 2), bty="n", las=1)
#' plot(o_mu1p2_NB, result="Prob", category=c(2, 1), bty="n", las=1)
#' 
#' }
#' @method plot ECFOCF
#' @export


# plot de la table ECF OCF ####

plot.ECFOCF <- function(x                  , 
                        ...                , 
                        result="CF"        , 
                        category=NA        , 
                        parameters="p1"    ,
                        legend=parameters  ,
                        period=1           , 
                        resultMCMC = NULL  , 
                        chain = "all"      , 
                        replicates = "all" ) {
  p3p <- list(...) # p3p=list()
  
  # result="CF"; category=NA; period=1; p3p=list(); x=NULL
  
  result <- tolower(result)
  result <- match.arg(result, choices = tolower(c("data", "dataOCF", 
                                                  "dataECF", 
                                                  "CF", 
                                                  "OCF", 
                                                  "ECF", 
                                                  "ECFOCF", 
                                                  "ECFOCF0", 
                                                  "probabilities", 
                                                  "period")))
  
  samples <- NULL
  if (!is.null(resultMCMC)) {
    if (chain == "all") {
      nbchains <- 1:(resultMCMC$parametersMCMC$n.chains)
      totalMCMC <- resultMCMC$resultMCMC[[1]]
      if (nbchains > 1) {
        for (j in 2:nbchains) {
          totalMCMC <- rbind(totalMCMC, resultMCMC$resultMCMC[[j]])
        }
      }
    } else {
      totalMCMC <- resultMCMC$resultMCMC[[chain]]
    }
    
    samples <- 1:(nrow(totalMCMC))
    if (replicates != "all") {
      replicates <- as.numeric(replicates)
      if (replicates <= max(samples)) {
        samples <- floor(seq(from=1, to=nrow(totalMCMC), length.out=replicates))
        # samples <- sample(x=samples, size=replicates, replace = FALSE)
      } else {
        samples <- sample(x=samples, size=replicates, replace = TRUE)
      }
    }
  } else {
    totalMCMC <- NULL
  }
  
  if (result=="data") {
    do.call(getFromNamespace("plot.TableECFOCF", ns="phenology"), 
            modifyList(list(x=x$data, 
                            period=period, 
                            cex.points=4, 
                            pch=19, 
                            col="black",
                            cex.axis=0.8, 
                            cex.labels=0.5, 
                            col.labels="red", 
                            show.labels=FALSE, 
                            show.0=FALSE, 
                            pch.0=4, 
                            cex.0=0.5, 
                            col.0="blue", 
                            show.scale = TRUE), p3p))
  }
  
  if (result=="ecfocf0") {
    if (all(is.na(category)) | (all(category == ""))) {
      do.call(getFromNamespace("plot.TableECFOCF", ns="phenology"), 
              modifyList(list(x=x$ECFOCF_0, 
                              period=period, 
                              cex.points=4, 
                              pch=19, 
                              col="black",
                              cex.axis=0.8, 
                              cex.labels=0.5, 
                              col.labels="red", 
                              show.labels=FALSE, 
                              show.0=FALSE, 
                              pch.0=4, 
                              cex.0=0.5, 
                              col.0="blue", 
                              show.scale = TRUE), p3p))
    } else {
      do.call(getFromNamespace("plot.TableECFOCF", ns="phenology"), 
              modifyList(list(x=x$ECFOCF_0_categories[[as.numeric(category)]], 
                              period=period, 
                              cex.points=4, 
                              pch=19, 
                              col="black",
                              cex.axis=0.8, 
                              cex.labels=0.5, 
                              col.labels="red", 
                              show.labels=FALSE, 
                              show.0=FALSE, 
                              pch.0=4, 
                              cex.0=0.5, 
                              col.0="blue", 
                              show.scale = TRUE), p3p))
    }
  }
  
  if (result=="ecfocf") {
    if (all(is.na(category)) | (all(category == ""))) {
      do.call(getFromNamespace("plot.TableECFOCF", ns="phenology"), 
              modifyList(list(x=x$ECFOCF, 
                              period=period, 
                              cex.points=4, 
                              pch=19, 
                              col="black",
                              cex.axis=0.8, 
                              cex.labels=0.5, 
                              col.labels="red", 
                              show.labels=FALSE, 
                              show.0=FALSE, 
                              pch.0=4, 
                              cex.0=0.5, 
                              col.0="blue", 
                              show.scale = TRUE), p3p))
    } else {
      do.call(getFromNamespace("plot.TableECFOCF", ns="phenology"), 
              modifyList(list(x=x$ECFOCF_categories[[as.numeric(category)]],
                              period=period, 
                              cex.points=4, 
                              pch=19, 
                              col="black",
                              cex.axis=0.8, 
                              cex.labels=0.5, 
                              col.labels="red", 
                              show.labels=FALSE, 
                              show.0=FALSE, 
                              pch.0=4, 
                              cex.0=0.5, 
                              col.0="blue", 
                              show.scale = TRUE), p3p))
    }
  }
  
  if (result=="cf") {
    # Si SE...
    if (!is.null(samples)) {
      CF_matrix <- universalmclapply(X=seq_along(samples), FUN = function(i) {
        par <- totalMCMC[samples[i], ]
        cf_e <- fitCF(x=par, fixed.parameters = x$fixed.parameters, data=x$data, 
                      itnmax = 0, hessian = FALSE)
        if (all(is.na(category)) | (all(category == ""))) {
          cf_e <- cf_e$CF
        } else {
          cf_e <- cf_e$CF_categories[[as.numeric(category)]]
        }
        return(cf_e)
      }, mc.cores = detectCores(), 
      clusterEvalQ=list(expr=expression(library(phenology))), 
      clusterExport=list(varlist=c("x", "totalMCMC", "category", "chain"), envir=environment()), 
      progressbar=TRUE)
      
      CF_matrix <- t(sapply(CF_matrix, FUN = function(x) x))
      
      q <- apply(X=CF_matrix, MARGIN=2, FUN=function(c) quantile(c, probs=c(0.025, 0.5, 0.975)))
      
      
      if (all(is.na(category)) | (all(category == ""))) {
        main="Clutch Frequency: All categories"
      } else {
        main=paste0("Clutch Frequency: Category ", as.character(category))
      }
      
      maxq <- ncol(q)
      
      do.call(plot_errbar, modifyList(list(x=1:maxq, 
                                           xlab="Clutch Frequency", 
                                           ylab="Density", 
                                           main=main, 
                                           y=q["50%", 1:maxq], 
                                           y.minus=q["2.5%", 1:maxq], 
                                           y.plus=q["97.5%", 1:maxq], 
                                           pch=19, 
                                           bty="n", 
                                           ylim=c(0, max(q["97.5%", 1:maxq])), 
                                           type="p", xaxt="n"), p3p)[c("x", "y", "type", "col", 
                                                                       "pch", "y.minus", "y.plus", "bty", 
                                                                       "main", "cex.axis", "bty", "las", 
                                                                       "xlab", "ylab", "xaxt", "xlim", "ylim")])
      
      
      axis(side = 1, at=1:maxq, cex.axis=unlist(modifyList(list(cex.axis=0.8), p3p)[c("cex.axis")]))
      
    } else {
      
      if (all(is.na(category)) | (all(category == ""))) {
        cf <- x$CF
        main="Clutch Frequency: All categories"
      } else {
        cf <- x$CF_categories[[as.numeric(category)]]
        main=paste0("Clutch Frequency: Category ", as.character(category))
      }
      
      do.call(plot, modifyList(list(x=1:length(cf), 
                                    xlab="Clutch Frequency", 
                                    ylab="Density", 
                                    main=main, 
                                    y=cf, type="h", xaxt="n"), p3p)[c("x", "y", "type", "col", 
                                                                      "main", "cex.axis", "bty", "las", 
                                                                      "xlab", "ylab", "xaxt", "xlim", "ylim")])
      axis(side = 1, at=1:length(cf), cex.axis=unlist(modifyList(list(cex.axis=0.8), p3p)[c("cex.axis")]))
    }
  }
  
  if ((result=="ocf") | (result=="dataocf")) {
    if (result=="ocf") {
      ylab="Density"
      if (all(is.na(category)) | (all(category == ""))) {
        ecfocf <- x$ECFOCF[, , period]
        main="Observed Clutch Frequency: All categories"
      } else {
        ecfocf <- x$ECFOCF_categories[[as.numeric(category)]][, , period]
        main=paste0("Observed Clutch Frequency: Category ", as.character(category))
      }
    } else {
      ylab="Frequency"
      ecfocf <- x$data[, , period]
      main="Observed OCF"
    }
    ocf <- rowSums(ecfocf, na.rm = TRUE)
    do.call(plot, modifyList(list(x=0:(length(ocf)-1), 
                                  xlab="Observed Clutch Frequency", 
                                  ylab=ylab, 
                                  main=main,
                                  y=ocf, type="h"), p3p)[c("x", "y", "type", "col", 
                                                           "main", "cex.axis", "bty", "las", 
                                                           "xlab", "ylab", "xlim", "ylim", 
                                                           "xaxt", "yaxt", "axes")])
  }
  
  if ((result=="ecf") | (result=="dataecf")) {
    if (result=="ecf") {
      ylab="Density"
      if (all(is.na(category)) | (all(category == ""))) {
        ecfocf <- x$ECFOCF[, , period]
        main="Estimated Clutch Frequency: All categories"
      } else {
        ecfocf <- x$ECFOCF_categories[[as.numeric(category)]][, , period]
        main=paste0("Estimated Clutch Frequency: Category ", as.character(category))
      }
    } else {
      ylab="Frequency"
      ecfocf <- x$data[, , period]
      main="Observed ECF"
    }
    ecf <- colSums(ecfocf, na.rm = TRUE)
    do.call(plot, modifyList(list(x=0:(length(ecf)-1), 
                                  xlab="Estimated Clutch Frequency", 
                                  ylab=ylab, 
                                  main=main, 
                                  y=ecf, type="h"), p3p)[c("x", "y", "type", "col", 
                                                           "main", "cex.axis", "bty", "las", 
                                                           "xlab", "ylab", "xlim", "ylim", 
                                                           "xaxt", "yaxt", "axes")])
    
  }
  
  if (result=="period") {
    if (all(is.na(category)) | (all(category == ""))) {
      y <- x$period
      main="All categories"
    } else {
      y <- x$period_categories[[as.numeric(category)]]
      main=paste("Category", category)
    }
    
    perr <- list(x=0:(length(y)-1), 
                 y=y, 
                 las=1, bty="n", 
                 ylab="Probability of nesting", xlab="Period", 
                 main=main)
    perr <- modifyList(perr, p3p)
    
    do.call(plot, perr)
  }
  
  if (result=="probabilities") {
    if (is.null(totalMCMC)) {
      
      namespar <- names(x$par)[grepl("^p|^a|^OTN", names(x$par))]
      if (identical(gsub("\\D", "", namespar), "")) {
        ncat <- 1
      } else {
        ncat <- max(as.numeric(gsub("\\D", "", namespar)), na.rm = TRUE)
      }
      
      par <- x$par[namespar]
      
      par_p <- invlogit(-par[grepl("^p", names(par))])
      if (length(par_p) == 1) if (ncat == 1) names(par_p) <- "p1" else {par_p <- rep(par_p, ncat); names(par_p) <- paste0("p", as.character(1:ncat))}
      par_a <- invlogit(-par[grepl("^a", names(par))])
      if (!identical(par_a, numeric(0))) {
        if (length(par_a) == 1) if (ncat == 1) names(par_a) <- "a1" else {par_a <- rep(par_a, ncat); names(par_a) <- paste0("a", as.character(1:ncat))}
      }
      par_OTN <- invlogit(-par[grepl("^OTN", names(par))])
      if (ncat == 1) names(par_OTN) <- "OTN1" else {
        if (length(par_OTN) == 1) {
          names(par_OTN) <- "OTN1"
        }
        par_OTN <- par_OTN[order(as.numeric(gsub("\\D", "", names(par_OTN))))]
        comp_par_OTN <- 1 - sum(par_OTN)
        names(comp_par_OTN) <- paste0("OTN", ncat)
        par_OTN <- c(par_OTN, comp_par_OTN)
      }
      
      total_columns <- length(par_p)+length(par_OTN)
      cn <- c(names(par_p), names(par_OTN))
      if (!identical(par_a, numeric(0))) {
        total_columns <- total_columns + length(par_a)*2
        cn <- c(cn, names(par_a), paste0("(a.p)", as.character(1:ncat)))
      }
      
      Total_par_prob <- matrix(data = NA, ncol=total_columns, nrow=1)
      colnames(Total_par_prob) <- cn
      
      par_p <- invlogit(-par[grepl("^p", names(par))])
      if (length(par_p) == 1) if (ncat == 1) names(par_p) <- "p1" else {par_p <- rep(par_p, ncat); names(par_p) <- paste0("p", as.character(1:ncat))}
      par_a <- invlogit(-par[grepl("^a", names(par))])
      if (!identical(par_a, numeric(0))) {
        if (length(par_a) == 1) if (ncat == 1) names(par_a) <- "a1" else {par_a <- rep(par_a, ncat); names(par_a) <- paste0("a", as.character(1:ncat))}
        par_ap <- par_a * par_p
        names(par_ap) <- paste0("(a.p)", as.character(1:ncat))
      }
      par_OTN <- invlogit(-par[grepl("^OTN", names(par))])
      if (ncat == 1) names(par_OTN) <- "OTN1" else {
        if (length(par_OTN) == 1) {
          names(par_OTN) <- "OTN1"
        }
        par_OTN <- par_OTN[order(as.numeric(gsub("\\D", "", names(par_OTN))))]
        comp_par_OTN <- 1 - sum(par_OTN)
        names(comp_par_OTN) <- paste0("OTN", ncat)
        par_OTN <- c(par_OTN, comp_par_OTN)
      }
      
      j <- 1
      if (any(grepl("^p", cn)))
        Total_par_prob[j, names(par_p)] <- par_p
      if (any(grepl("^OTN", cn)))
        Total_par_prob[j, names(par_OTN)] <- par_OTN
      if (!identical(par_a, numeric(0))) {
        Total_par_prob[j, names(par_a)] <- par_a
        Total_par_prob[j, paste0("(a.p)", as.character(1:ncat))] <- par_ap[paste0("(a.p)", as.character(1:ncat))]
      }
      main_ec <- "No uncertainty is shown"
    } else {
      # Dans parameters, j'ai les paramètres à montrer
      
      namespar <- colnames(totalMCMC)[grepl("^p|^a|^OTN", colnames(totalMCMC))]
      if (identical(gsub("\\D", "", namespar), "")) {
        ncat <- 1
      } else {
        ncat <- max(as.numeric(gsub("\\D", "", namespar)), na.rm = TRUE)
      }
      
      # J'utilise totalMCMC[samples[i], ]
      par <- totalMCMC[samples[1], namespar]
      
      par_p <- invlogit(-par[grepl("^p", names(par))])
      if (length(par_p) == 1) if (ncat == 1) names(par_p) <- "p1" else {par_p <- rep(par_p, ncat); names(par_p) <- paste0("p", as.character(1:ncat))}
      par_a <- invlogit(-par[grepl("^a", names(par))])
      if (!identical(par_a, numeric(0))) {
        if (length(par_a) == 1) if (ncat == 1) names(par_a) <- "a1" else {par_a <- rep(par_a, ncat); names(par_a) <- paste0("a", as.character(1:ncat))}
      }
      par_OTN <- invlogit(-par[grepl("^OTN", names(par))])
      if (ncat == 1) names(par_OTN) <- "OTN1" else {
        if (length(par_OTN) == 1) {
          names(par_OTN) <- "OTN1"
        }
        par_OTN <- par_OTN[order(as.numeric(gsub("\\D", "", names(par_OTN))))]
        comp_par_OTN <- 1 - sum(par_OTN)
        names(comp_par_OTN) <- paste0("OTN", ncat)
        par_OTN <- c(par_OTN, comp_par_OTN)
      }
      
      total_columns <- length(par_p)+length(par_OTN)
      cn <- c(names(par_p), names(par_OTN))
      if (!identical(par_a, numeric(0))) {
        total_columns <- total_columns + length(par_a)*2
        cn <- c(cn, names(par_a), paste0("(a.p)", as.character(1:ncat)))
      }
      
      Total_par_prob <- matrix(data = NA, ncol=total_columns, nrow=length(samples))
      colnames(Total_par_prob) <- cn
      
      
      for (j in samples) {
        par <- totalMCMC[samples[j], namespar]
        
        par_p <- invlogit(-par[grepl("^p", names(par))])
        if (length(par_p) == 1) if (ncat == 1) names(par_p) <- "p1" else {par_p <- rep(par_p, ncat); names(par_p) <- paste0("p", as.character(1:ncat))}
        par_a <- invlogit(-par[grepl("^a", names(par))])
        if (!identical(par_a, numeric(0))) {
          if (length(par_a) == 1) if (ncat == 1) names(par_a) <- "a1" else {par_a <- rep(par_a, ncat); names(par_a) <- paste0("a", as.character(1:ncat))}
          par_ap <- par_a * par_p
          names(par_ap) <- paste0("(a.p)", as.character(1:ncat))
        }
        par_OTN <- invlogit(-par[grepl("^OTN", names(par))])
        if (ncat == 1) names(par_OTN) <- "OTN1" else {
          if (length(par_OTN) == 1) {
            names(par_OTN) <- "OTN1"
          }
          par_OTN <- par_OTN[order(as.numeric(gsub("\\D", "", names(par_OTN))))]
          comp_par_OTN <- 1 - sum(par_OTN)
          names(comp_par_OTN) <- paste0("OTN", ncat)
          par_OTN <- c(par_OTN, comp_par_OTN)
        }
        
        
        if (any(grepl("^p", cn)))
          Total_par_prob[j, names(par_p)] <- par_p
        if (any(grepl("^OTN", cn)))
          Total_par_prob[j, names(par_OTN)] <- par_OTN
        if (!identical(par_a, numeric(0))) {
          Total_par_prob[j, names(par_a)] <- par_a
          Total_par_prob[j, paste0("(a.p)", as.character(1:ncat))] <- par_ap[paste0("(a.p)", as.character(1:ncat))]
        }
      }
      main_ec <- "95% quantiles using Bayesian MCMC"
    }
    qprob <- apply(Total_par_prob, MARGIN = 2, FUN = function(x) {quantile(x, probs = c(0.025, 0.5, 0.975))})
    qprob <- qprob[, parameters, drop=FALSE]
    
    perr <- list(x=1:ncol(qprob),  
                 y=qprob["50%", ], 
                 y.minus=qprob["2.5%", ], 
                 y.plus = qprob["97.5%", ], 
                 xlim=c(0.5, ncol(qprob)+0.5), 
                 ylim=c(0,1), 
                 las=1, bty="n", xaxt="n", 
                 ylab="Probability of capture", xlab="Categories", 
                 main=main_ec)
    perr <- modifyList(perr, p3p)
    
    do.call(plot_errbar, perr)
    
    if (is.null(p3p$xaxt)) p3p$xaxt <- "r"
    
    if (p3p$xaxt != "n") {
      segments(x0=1:ncol(qprob), 
               y0=-0.01, y1=-0.04, xpd=TRUE)
      cex <- p3p[["cex.axis"]]
      y <- p3p[["y.axis"]]
      if (is.null(cex)) cex <- 1
      if (is.null(y)) y <- -0.08
      do.call(text, modifyList(list(x = 1:ncol(qprob), 
                                    y=y, 
                                    cex=cex, 
                                    labels = legend, 
                                    xpd=TRUE), p3p[c("srt", "labels")]))
    }
  }
}